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ABSTRACT 

We present a global linear stability analysis of nuclear fuel accumulating on 
the surface of an accreting neutron star and we identify the conditions under 
which thermonuclear bursts are triggered. The analysis reproduces all the recog- 
nized regimes of hydrogen and helium bursts, and in addition shows that at high 
accretion rates, near the limit of stable burning, there is a regime of "delayed 
mixed bursts" which is distinct from the more usual "prompt mixed bursts." 
In delayed mixed bursts, a large fraction of the fuel is burned stably before the 
burst is triggered. Bursts thus have longer recurrence times, but at the same time 
have somewhat smaller fluences. Therefore, the parameter a, which measures the 
ratio of the energy released via accretion to that generated through nuclear reac- 
tions in the burst, is up to an order of magnitude larger than for prompt bursts. 
This increase in a near the threshold of stable burning has been seen in obser- 
vations. We explore a wide range of mass accretion rates, neutron star radii and 
core temperatures, and calculate a variety of burst properties. From a prelim- 
inary comparison with data, we suggest that bursting neutron stars may have 
hot cores, with T core > 10 K, consistent with interior cooling via the modified 
URCA or similar low-efficiency process, rather than T core ~ 10 7 K, as expected 
for the direct URCA process. There is also an indication that neutron star radii 
are somewhat small < 10 km. Both of these conclusions need to be confirmed by 
comparing more careful calculations with better data. 
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1. Introduction 

When gas accretes onto a neutron star, nuclear reactions often occur in an unstable 
fashion (Hansen & van Horn 1975), leading to thermonuclear explosions which are called 
Type I X-ray bursts. These bursts were first observed by Grindlay et al. (1976), and have 
been studied intensively for many years (see van Paradijs et al. 1988; Lewin et al. 1993; 
Strohmayer et al. 1998; Cornelisse et al. 2003, for summaries of the observations). Recently, 
there has been renewed excitement in the field, following the discovery of high frequency 
oscillations in burst light curves (e.g., Strohmayer et al. 1996; Strohmayer 2001a; van Straaten 
et al. 2001; Muno et al. 2001). 

The physics of Type I bursts has been widely studied, and the broad features of the 
phenomenon are understood theoretically (Woosley & Taam 1976; Joss 1977; Taam & Pick- 
lum 1978; Paczynski 1983a; Fujimoto et al. 1981; Fushiki & Lamb 1987a; Taam et al. 1996; 
Bildsten 1998). Models show that there are three kinds of bursts, depending on whether 
hydrogen or helium burning dominates. At high mass accretion rates M, hydrogen is only 
partially burned before a burst is triggered by unstable helium burning. The result is a mixed 
burst in which both hydrogen and helium burn explosively. At somewhat lower values of M, 
all the hydrogen is consumed before the helium instability is triggered. This leads to a pure 
helium burst. For yet lower M, hydrogen itself burns unstably, giving a hydrogen burst. The 
ranges of M corresponding to the different regimes have been worked out approximately. 

Two very different approaches have been pursued for theoretically modeling the burst 
phenomenon. In one approach, one simulates the physics of the accreting gas with a fully 
time-dependent code that includes a large network of nuclear reactions and sophisticated 
thermodynamics (e.g., Joss 1978; Taam & Picklum 1979; Joss & Li 1980; Taam 1982; Fu- 
jimoto et al. 1987b; Taam et al. 1993). Such studies are essential for understanding the 
details of the thermonuclear explosion; indeed, the one-dimensional simulations of the past 
have now been generalized to two-dimensional and even three-dimensional simulations (the 
FLASH effort at Chicago, Zingale et al. 2001) which follow the physics of bursts in exquisite 
detail. Such simulations are, however, not very convenient for parameter surveys or for 
detailed comparisons of theoretical predictions with observational burst statistics. 

An alternate approach, where one focuses on the thermonuclear instability that triggers 
the burst rather than on the burst itself, has been popular (e.g., Hansen & van Horn 1975; 
Ergma & Tutukov 1980; Fujimoto et al. 1981; Taam 1982; Paczynski 1983a; Fushiki & Lamb 
1987a; Fujimoto et al. 1987a; Cumming & Bildsten 2000). Here one treats the accumulating 
layer of gas on the surface of the neutron star as a quasi-equilibrium system whose properties 
vary slowly with time. One first solves for the equilibrium structure and then analyses the 
stability of the gas layer by considering the effect of small perturbations on the underlying 
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equilibrium. If the perturbations grow with time, one says that the system is unstable, 
presumably resulting in a Type 1 X-ray burst. 

The second approach, though not as general as the first, allows one to explore a large 
range of parameter space and to study how the burst phenomenon is affected by variations 
in control parameters such as the mass accretion rate, the surface gravity of the star, etc. 
However, the analyses that have been published so far in the literature employ fairly simple- 
minded criteria to decide exactly when an accreted layer becomes unstable and are, therefore, 
not very accurate or complete. The motivation of the present paper is to develop a more 
rigorous stability analysis for accreting neutron stars. Our hope is (i) to verify the results of 
the earlier methods, (ii) to explore additional phenomena that may have been missed previ- 
ously, and (iii) to put the theory on a rigorous footing to enable quantitative comparisons 
with observations. We reported some early results of this work in Narayan & Heyl (2002). 

We begin the paper in §2 with a description of our model of the accreted layer and 
the method we employ to calculate equilibria and to study their stability. We follow this 
in §3 with an exploration of sequences of equilibria, and in §4 with a discussion of the 
stability properties of the equilibria. Through this analysis, we reproduce the three previously 
recognized regimes of burst activity, namely mixed bursts, helium bursts, and hydrogen 
bursts. However, we find that mixed bursts themselves come in two kinds: prompt mixed 
bursts and delayed mixed bursts. The latter category has not been recognized previously. In 
§5 we present results for neutron stars of various radii and core temperatures, and explore 
a wide range of mass accretion rates. We calculate a number of burst observables such 
as the recurrence time, the burst duration, and the dimensionless parameter a (Eq. 35) 
which is widely used in the burst literature. In §6, we compare our theoretical formalism 
to methods previously published in the literature. We also compare our predictions with 
selected observations. From the latter, we obtain preliminary results on the likely core 
temperatures and radii of bursting neutron stars. We conclude in §7 with a summary. 



2. The Model 

2.1. Governing Equations 

We assume that gas accretes on the surface of a compact spherical star of mass M 
and radius it! at a rate E (mass per unit area per unit time). We consider all quantities 
to be functions of S, the column density (mass per unit area) measured from the top of 
the accreted layer. We use partial derivatives d/dt and <9/<9X to signify Eulerian time and 
"spatial" derivatives at a fixed S, and d/dt to represent the "Lagrangian" time derivative 
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following a parcel of accreted gas (see Fushiki & Lamb 1987a): 

d _ d ■ d 

We consider only hydrogen and helium burning, and so we describe the composition of the 
gas in terms of the hydrogen mass-fraction X, the helium mass fraction Y, and the heavy 
element fraction Z = 1 — X — Y . These quantities start off with values X out , Y out , Z ont at 
£ = 0, corresponding to the composition of the gas initially falling on the neutron star, and 
evolve with increasing depth as a result of nuclear burning. Because if -burning is mostly 
done via the CNO-cycle, it is necessary to know what fraction of Z is in CNO elements (see 
Eq. 26). For this, we assume that about 80% of the initial Z ont is in CNO, as appropriate for 
solar composition (Allen 2000), and that all the Z produced via helium burning is entirely 
in CNO; thus, we take 

Zcno = 0.8Z out + (Z — Z out ). (2) 

The time evolution of the accreting gas is described by a set of five coupled partial 
differential equations: 



dP 
d£ 

dT 3kF 

d£ 16aT 3 ' 
OF ds 

5E = - T ^-( e H + e He ), (5) 



9, (3) 

(4) 



dX _en_ 
~dW ~ 

dF e H 6hc 



dt E* Hc 



(6) 
(7) 



Table 1 gives the definitions of the various symbols. We assume that the accreted layer 
is thin relative to the stellar radius, and so we take the gravitational acceleration g to be 
independent of E: 

. GM , / 2GMY 1 ' 2 

where z is the gravitational redshift. Note that g and all other quantities in equations (3-7) 
are measured in the local frame of the gas. In this spirit, E is the baryonic mass added per 
unit surface area per unit local time. 

Equation (3) describes the condition of hydrostatic equilibrium. By making use of this 
equation rather than the full momentum equation, we filter out sound waves and focus on 
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Table 1. Definition of symbols: 



Symbol 



Meaning 



R, M, R s 

9 z 
t 

Lgjco -^Eddj ^acc 

s 

y 

^ layer 

^ layer, crit 

^diff) S max 

t, P, p, T, s 
F 

X, Y,Z 

-ZcNO 

Po, Tq] pi, Ti, 
Tout; ^out? etc. 
F 

± nuc 

-^out) /nuc 
± acc 

TJayer; -^layer; etc. 

T T 

- 1 max ; - 1 core 

Tdiff 

K, K 

CH,He 

Z?* 
H,He 

7 

7aco ^acc 
^rcc 

Eft, i?He 
^H+He, ^He 

a 



stellar radius, stellar mass, Schwarzschild radius: 2GM/c 2 
gravitational acceleration, redshift at the surface 
mass accretion rate per unit area 

accretion luminosity, Eddington luminosity, Z acc = L acc /L-Edd 

surface mass density measured from the surface; independent variable in eqs. 

E of the accreted layer 

critical Si aycr at which a burst is triggered 

Maximum integration depth for the equilibrium, and for perturbations 

time, pressure, mass density, temperature, entropy per unit mass 

outbound energy flux 

mass fractions of hydrogen, helium, metals 

mass fraction of CNO nuclei 

values in equilibrium; perturbations 

values at the surface 

expected surface flux if entire fuel is steadily burned 

escaping flux due to nuclear reactions and compression, / out = F ont /F nnc 

persistent flux from the surface due to accretion: H(?zj(\ + z) 

values at the bottom of the accreted layer 

maximum T in layer, core temperature 

thermal diffusion time 

opacity, conductivity 

energy-generation rate per unit mass for H, He burning 
total nuclear energy released per unit mass of H, He burned 
growth rate of mode, 7 = 9^(7) + 2^(7) 

7 acc = S/Ei ayer , accretion rate; t acc = S laycr /S, accretion time 
burst recurrence time: (1 + z)£i ayeriCrit /£ 
fluence in burst from burning hydrogen, helium 
effective burst duration: (_Eh + E^ e ) / L-^dd, E^ e /L-Edd 
ratio of accretion energy to nuclear energy: (t rC c/^H+He)Lcc 
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variations that occur on a much longer time scale than the sound-crossing time. Equation (4) 
describes energy transfer by radiation and conduction. Equation (5) is the energy conserva- 
tion equation, and equations (6) and (7) describe the evolution of the hydrogen and helium 
mass-fractions as a result of nuclear burning. Equations (3)-(7) require expressions for the 
pressure P, the opacity k, the entropy s, and the energy-generation rates for hydrogen and 
helium-burning, en, chg- These are discussed in §2.3. 

The calculations proceed as follows. We start with a bare neutron star and follow the 
properties of the accreted layer as gas piles up. For a given column density Si ayer of the 
accreted layer, we calculate the following. First, we solve for the quasi-equilibrium state of 
the system. To do this, we set the Eulerian time derivative d/dt to zero in equations (3-7) 
and consider the following set of ordinary differential equations (see Fushiki & Lamb 1987a) 



dp 

dE 


= 9, 


(9) 


dT 
dE 


3kF 
16aT 3 ' 


(10) 


dF 
dE 


= -ST— - (e H + e He ), 


(11) 


dx 

dE 




(12) 


E— 
dE 


CH ^He 
TP* T?* 


(13) 



We solve these equations with boundary conditions (described in §2.2) to obtain the run 
of density Po(^), temperature T (E), etc., in quasi-steady state. Note that, because Si aycr 
increases steadily with time, the layer at any given time is not strictly in equilibrium, and 
the above steady state equations are not precisely valid. However, since Ei ayer increases only 
slowly with time (on the accretion time scale defined in eq. 17 below), and since many of 
the physical processes in the layer have shorter characteristic time scales, this is a reasonable 
approximation. 

Having calculated the steady state equilibrium solution, we carry out a linear pertur- 
bation analysis. This is the principal contribution of our work, and represents a significant 
advance over previous studies. For the perturbation analysis, we assume that the various 
physical quantities are functions of E and t in the form 

p(E,f) = po(S) + exp(7t) Pl (E), (14) 
T(E, t) = T (E)+exp( 7 t)T 1 (E), etc., (15) 

where po, T represent the solutions obtained from solving the steady state equations de- 
scribed above, and p\ <C p , T\ <C T , are small linear perturbations. The frequency 7 



-7- 



represents the growth rate of the perturbations. We substitute the perturbed solution into 
the original time-dependent equations (3-7) and linearize in the usual way to obtain a set 
of ordinary differential equations (in S) for the first-order quantities pi, T 1; etc. These 
equations are written down in Appendix A and some of their properties are discussed there. 
We solve the linearized perturbation equations with appropriate boundary conditions and 
thereby obtain 7, which plays the role of an eigenvalue. 

In general, there are many solutions for 7 for a given steady state solution; some values 
of 7 are real and some are complex. If any solution for 7 has a real part sufficiently large 
compared to the accretion rate — see equation (21) below for a precise statement of what 
the criterion is — then we say that the layer is unstable. In this case, we identify 9^(7) with 
the growth rate of the instability in the system. If no solution for 7 satisfies equation (21), 
then we consider the system to be stable. Because 7 is in general complex, all the quantities 
in the perturbation equations are complex and must be handled with complex arithmetic (in 
contrast to the steady state equations which involve purely real quantities). 

We carry out the above two stages of calculations, namely steady state and linear 
perturbation analysis, for each S layer as matter accumulates on the star. If no instability is 
found for any choice of Si ay cr up to a very large value, typically 10 13 — 10 14 g cm~ 3 — see 
Taam & Picklum (1978) and Brown & Bildsten (1998) for a discussion of carbon flashes which 
occur at yet larger column densities — then we say that the system is stable to bursts. If, for 
some value of Si ayer , we do obtain an instability, then we say that the system will undergo a 
burst when it accumulates this much gas on its surface. 

2.2. Boundary Conditions 

The solution of the five differential equations (9-13) for the steady state requires five 
boundary conditions. Four are applied on the outside, at the photosphere (where the optical 
depth is taken to be 2/3), and one is applied below the accreted layer, deep inside the star. 

The outer boundary conditions are very similar to the ones employed in Narayan & 
Heyl (2002). The surface values of X and Y are set equal to the corresponding values of the 
accreting gas, namely X = X out , Y = Y out . We have used a solar composition, X out = 0.7, 
^out — 0.28 (Allen 2000, Table 3.1), in all the calculations reported in this paper. Next, a 
particular value is assumed for the escaping flux F out from the accreted layer. This flux is 
the result of nuclear burning and compression, and its precise value is determined only after 
applying the inner boundary condition, as explained below. The surface temperature T out is 
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then obtained by the condition 

aT 4 ut = F out + F acc , F acc = ±c 2 z{l + z), (16) 

where F acc is the gravitational energy released by the accreting gas, and we have included 
the appropriate gravitational redshift factor so that all quantities are calculated in the local 
frame. The above relation for T out is approximate, but we have confirmed that the results 
are very insensitive to the precise choice of T out . The fourth boundary condition is obtained 
from the radiative transfer equation. Given an assumed value for F out and the condition 
P = at E = 0, this equation directly gives the density profile p(E). 

As described above, the solution at the surface is completely specified once a value of 
F out is assumed. The unknown value of -F out is determined by requiring the solution to satisfy 
an inner boundary condition. For this, we assume that the accreting star has a specified 
core temperature T core and we require the solution of the steady state equations to match 
this temperature. The key issue is where exactly to do the matching. Our approach is as 
as follows. Associated with a given accreted column of depth £i ayer , there is a characteristic 
accretion time 

tacc = (17) 

We integrate the steady state equations from the surface down to the bottom of the accreted 
layer at E layer , and then we integrate further into the stellar substrate to a depth £ diff such 
that the diffusion time from £i aycr to S diff is equal to twice t acc (the factor of 2 is arbitrary 
and was selected after some numerical experiments): 

7diff(£diff) — 2t a cc- (18) 

For applying this condition, we need an estimate of the thermal diffusion time Tdiff(E) for any 
choice of £ inside the star. We obtain this by integrating a separate differential equation: 

dr diff = 9/tA: B T(E - E laycr ) 

dE ~ l28aT*fxm u ' [ ' 

The form of this equation and the numerical coefficient (which is based on a simple toy model) 
are obtained by treating the energy equation (11) as a diffusion problem. The precise details 
are unimportant for the final results. 

For the calculations presented in this paper, we pick "reasonable" values for T corc , trying 
a range that is likely to bracket the true value. In §5.5 we will estimate the core temperature 
as a function of accretion rate for two kinds of neutrino cooling in the core. If we include both 
the inward directed flux at the bottom of the layer that we calculate and the energy from 
nuclear reactions in the deep crust (Brown et al. 1998), we find that the latter contribution 
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dominates. This yields a simple relationship between the accretion luminosity and the core 
temperature. 

The above treatment of the inner boundary condition is superior to the method used in 
Narayan & Heyl (2002, and most other previous studies, see §6.1), where the temperature 
at the bottom of the accreted layer was set equal to T corc . That is, in that calculation, 
the substrate was assumed to be isothermal immediately below the accreted layer. By 
integrating down to a couple of diffusion lengths into the substrate, we believe our present 
approach is physically better motivated. However, even with this approach, we are effectively 
assuming that the core below Edig is perfectly isothermal — hence the use of the term "core 
temperature" — which is again an approximation. For a system that has been accreting and 
bursting for a long time, there is expected to be a small time-averaged flux flowing into the 
core even inside of £ = Edis- In addition, there is a larger flux from deep crustal reactions 
(Brown et al. 1998). Both of these fluxes will induce a temperature gradient, so that the 
temperature in the crust, which is relevant for bursts, will not be equal to the temperature 
T core deep inside the neutron star (set for instance by neutrino cooling, see Brown 2000 for 
a detailed analysis). This effect may not be very serious, since the matter inside Ediff is 
highly degenerate and very conductive. Nevertheless, we mention the point here because the 
effect of the approximation is presently not fully understood. We believe that the results 
we present in this paper for T corc ~ 10 8 K will be hardly affected because the burning layer 
itself has a temperature of this order. However, for very cold cores, e.g., T core ~ 10 7 K, the 
approximation may have a more serious effect. 

The above discussion pertains to the equilibrium solution. The linear perturbation 
equations have similar boundary conditions. At the surface, the first-order perturbations 
X 1 and Y 1 vanish, and we assume an arbitrary value for the flux perturbation the latter 
choice serves as the overall normalization of the perturbed eigenfunction. From the perturbed 
flux, the corresponding temperature perturbation 7\ is readily obtained via equation (16), 
and finally the density perturbation p\ is calculated from the radiative transfer equation. 

At the bottom, we set the temperature perturbation Ti(E max ) at a prescribed depth 
E max equal to zero. In analogy with what we did for the steady state solution, we determine 
E max by the condition that the diffusion time down to this depth should be equal to twice 
the mode time scale t mo d e : 

2 

T difr(E max ) = 2t mode = -, — r, (20) 

ItI 

where |7| is the modulus of the complex eigenvalue 7. We only consider modes that grow 
faster than the accretion time, see below; therefore, £ max is always smaller than Ediff- The 
condition Ti(£max) = provides the final boundary condition that enables us to solve for 
the eigenmode and the eigenvalue 7. 
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As already mentioned, there are many solutions for 7. We concentrate on modes that 
grow fast enough to be interesting, specifically modes that satisfy 

1 S 

3^(7) > 0mode7acc, fi'mode = 3, 7 acc = = — . (21) 

f-acc ^layer 

The factor 3 is arbitrary, but reasonable. For an instability to have any noticeable effect on 
the system, it needs to grow in a time shorter than the lifetime t acc of the system. Also, 
our approach of treating the accreted layer as a quasi-steady system, and analysing the 
perturbations as if they occur in a time-steady system, is valid only if the time scale of the 
perturbations is sufficiently small compared to t acc . For both reasons, the choice, g mo de = 3, 
in equation (21) seems appropriate. 



2.3. Auxiliary Prescriptions 

The solution of the set of differential equations described in §2.1 requires knowledge of 
the thermodynamic and other physical properties of the gas. We describe here the particular 
prescriptions we have used for the calculations presented in this paper. 



2.3.1. Equation of State 

The pressure is assumed to be supplied by photons, nuclei and electrons. For the photons 
we use the blackbody formula and for the nuclei we assume an ideal non-degenerate gas. In 
the case of the electrons, we write the pressure as the quadrature sum of two terms, one 
equal to the pressure of an ideal non-degenerate gas and the other equal to the pressure of 
a zero-temperature degenerate electron gas. Thus we take (see Paczynski 1983b) 

P = + P„ d + [P c 2 nd + P c 2 d ] 1/2 , (22) 

_ pk B T _ pk B T 

- r nuc,nd j r e,nd i V Z " J / 

where the molecular weights /i nuc and p e are determined in the standard way (e.g. Clayton 
1983). In determining /i nuc , A*e, we assume that a fraction 0.2Z ont of the heavy elements in the 
accreted layer consists of 56 Fe and the rest (what we have called -Zcno, see eq. 2) consists 
of 14 N (as a surrogate for CNO elements). The latter choice is an approximation. The 
initial CNO elements in the accreted gas consist mostly of 12 C and 16 0. During H-burning 
via the CNO cycle, the composition is mostly 14 and 15 0, while the composition reverts 
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back to 12 C and 16 after the CNO cycle is done. Considering the other simplifications we 
have employed, approximating the composition as 14 N for the purposes of calculating the 
molecular weights seems reasonable. We take the stellar substrate below the accreted layer 
to be made of pure 56 Fe. 

For the zero-temperature degenerate electron pressure Pdcg, we use an exact expression 
that is valid for all Fermi energies (Shapiro & Teukolsky 1983). We do not include Coulomb 
corrections on the pressure. Separating the electron pressure into a non-degenerate part and 
a zero-temperature degenerate part and taking the quadrature sum is, again, an approxi- 
mation. The separation works well in various asymptotic limits and is reasonably accurate 
even in the transition regime between non-degeneracy and degeneracy (Paczynski 1983b). 
We believe that the approximation is adequate for our purposes since the transition zone is 
usually fairly narrow in S. 

The entropy term in equation (5) is important since it is the origin of the compressional 
flux. We write the entropy as the sum of contributions from each nuclear species and the 
electrons. The entropy per unit mass for species % takes the form 



Nik 



,/, t n 3,_ 5 3 , / 2nrriik 
ln ( pj V,) + _ ln T + - + -, n (^ 



(24) 



where JVj is the number of particles of the particular species per unit mass of the gas, rrii 
is the mass of each particle, and h is Planck's constant. This expression corresponds to an 
ideal gas. In the case of the electrons, we use the above entropy so long as the quantity is 
positive, and replace it with zero when the expression becomes negative. This allows us to 
handle both the non-degenerate and degenerate limits adequately. 



2.3.2. Opacity 

We include radiative and conductive energy transfer, and model the opacity k as 

111 

- = + • (25) 

^ ^rad ^cond 

For the radiative opacity « ra d, we use the formulae given in Appendix A of Schatz et al. 
(1999), using the analytic formulae of Antia (1993) to calculate the electron chemical poten- 
tial. 

To calculate the conductive opacity we have used the software of Potekhin (1999). Up 
to a density of 10 9 g cm -3 , we include the effects of impurities among the nuclei. Above this 
density we assume that the material is pure, and we interpolate between the two regimes. 
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The density of the fuel layer is nearly always less than 10 9 g cm 3 , except for some rare 
helium bursts. 



2.3.3. Nuclear Energy Generation Rates 

For e H , we include the pp chain and the CNO cycle. Because the temperature of the 
burning material is on the order of 10 8 K, the CNO cycle usually dominates. We include 
fast-CNO burning, saturated CNO burning and electron capture reactions, as described in 
Mathews & Dietrich (1984) and Bildsten & Cumming (1998), and write the energy-generation 
rate as 

eH,CNO = 4£HrcNO— ^j— , (26) 

where tcno is the rate of reactions per CNO nucleus, and we have assumed that in equilibrium 
the majority of the CNO nuclei are 14 N (Clayton 1983). For the reaction rate, we assume 



/ 1 

?~CNO = h 



\Ti3 862.0 s 



1 



+ 



7i3/ \ r i3 + r i4 + 278.2 s / V 862 - s / V r i4 + 1038.0 s 

(27) 

Here, r 13 is the lifetime of 13 N against the reaction 13 N(p, 7) 14 (Mathews & Dietrich 1984), 

r 13 = (Xp/gcnT 3 )' 1 3.35 x 10 7 T~ 2/3 exp (-15.202 T~ 1/3 - 0.8702 T 9 2 ) 

x (l + 0.027 T 9 1/3 + 0.9 T 9 2/3 + 0.173 T 9 + 4.61 T 9 4/3 + 2.26 T 9 /3 ) (28) 
+ 3.03 x 10 5 T 9 _3/2 exp (-6.348 T 9 l ) 



-l 

s, 



where T 9 = T/(10 9 K), and r 14 is the lifetime of 14 N against the reaction 14 N(p, 7) 15 (e.g. 
Paczynski 1983a), 

7i4 = 3.1 x 10 10 (Xp/gcm" 3 )~ 1 T 6 2/3 exp (l52.313 T~ 1/3 ) s, (29) 

where T 6 = T/(10 6 K). The time scale 862.0 s refers to the beta-decay timescale of 13 N, 
278.2 s is the sum of the beta-decay timescales of 14 and 15 0, and 1038.0 s is the sum of 
the beta-decay timescales of 13 N and 15 0. In deriving the above rates we have assumed that 
the various species have reached their equilibrium abundances. 

Figure 1 shows the variation of en and ch c with temperature for some typical densities. 
For low and moderate temperatures, the hydrogen-burning rate is a steeply increasing func- 
tion of temperature, and it is this steep dependence that drives a thermonuclear instability. 
For a temperature greater than about 10 7 ' 8 — 10 7 ' 9 K, however, hydrogen-burning switches 
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rather abruptly to the saturated burning regime. Beyond this temperature, hydrogen- 
burning is stable. As we discuss in §3, this change has a noticeable effect on the sequence of 
equilibria. We have corrected the proton capture rates for screening using the formulae of 
Dewitt et al. (1973) for the non- resonant reactions and the formulae of Itoh et al. (2003) for 
the resonant reaction ( 13 N(p, 7) 14 0). Screening increases the reaction rates at temperatures 
where the CNO cycle is not saturated. 

For eHe, we use the fitting formula given in equations (4.7), (4.8a) and (4.8b) of Fushiki 
& Lamb (1987b) which include screening. We introduce a smooth transition between the 
various regimes defined by these authors in order to more faithfully reproduce the numerical 
results they have depicted in their Fig. 3. Our fitting results are shown in Fig. 1. 



3. Equilibrium Solutions 

Paczynski (1983a) has presented a very helpful analysis of the stability of nuclear burning 
on the surface of a compact star. His model involves numerous simplifications: he uses a one- 
zone approximation, he considers only helium burning, and he assumes an inner boundary 
condition on the flux rather than on the temperature. Nevertheless, many of the insights 
he has obtained via his simple analysis carry over to our more detailed work. In particular, 
following his work, we have found it very helpful to consider equilibria in the -F out -E layer 
plane. Appendix A discusses why the insights from Paczyhski's analysis apply to our more 
complicated model. 

Figures 2 and 3 show sequences of equilibria of the accreted layer for a 1.4M Q neutron 
star with a radius of 10.4 km and a core temperature of 10 8 K. The different panels correspond 
to different mass accretion rates M, parameterized by the accretion luminosity L acc measured 
at infinity, 

L acc = Mc 2 ^ — ^ = /acc-^Edd, (30) 

where L Edd = ^GMcj n es , with k, cs = 0.4 cm 2 g _1 , is the Eddington luminosity measured at 
infinity. The local surface mass accretion rate E is related to M by 

M 

~ ' (1 + *). (31) 



4ttR 2 



The mass accretion rate M, or equivalently E, is a key parameter that determines the nature 
of bursts. Since it is however not directly measured, we prefer to give all our results in terms 
of the dimensionless luminosity l acc = L acc /LEdd- In doing this, we assume that the accretion 
is radiatively efficient and satisfies the relation given in equation (30). 
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Fig. 1. — Nuclear reaction rates for hydrogen and helium, plotted as a function 
of temperature. From bottom to top, the four curves correspond to densities of 
10 5 , 10 5 ' 5 , 10 6 , 10 6 ' 5 gem -3 . The hydrogen-burning curve is dominated by the pp chain 
at low temperatures, the CNO cycle at intermediate temperatures, and saturates above 
about 10 7 - 8 — 10 7 - 9 K. The dashed lines trace the results if one ignores screening. The helium 
rates are taken from Fushiki & Lamb (1987b). 
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Fig. 2. — The points trace the locus of values of the surface density of the accreted layer 
Si aycr and the normalized escaping flux / out = F out /F nuc which satisfy the outer and inner 
boundary conditions. The calculations assume a neutron star of mass M = 1.4 M Q , radius 
R = 10°' 4 i? s = 10.4 km, and core temperature T core = 10 s K. The four panels correspond to 
four relatively large accretion luminosities. An open circle indicates that the corresponding 
equilibrium is stable. A filled circle indicates that the equilibrium is unstable with a real 
mode frequency 7, while a star indicates an unstable equilibrium with complex 7, i.e., an 
overstability. 
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Fig. 3. — Similar to Fig. 2, but for four lower accretion luminosities. 
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The eight panels in Figs. 2, 3 correspond to accretion luminosities log(7 acc ) = —0.5, 
—0.75, — 1, —1.25, —1.5, —1.75, —2.75 and —3 in Eddington units. In each panel, the 
horizontal axis shows the escaping flux / out from the accreted layer, normalized by the 
maximum nuclear burning energy available in the accreting gas: 

/out = -jj, — , F nuc = £ [X out Eft + (X out + Y out )Eft e ] . (32) 

-*nuc 

The vertical axis shows the column density of the layer, £i ayer - 

The most obvious feature of the various panels is that the equilibria do not form mono- 
tonic sequences in the F out -Y,\ aycr plane. At a given S laycr , there can be one, three, or even 
five, distinct solutions for F out . (There are no cases of five solutions in the sequences shown 
here, but it is fairly common when the core temperature is lower, e.g., 10 75 K). Because the 
equations are nonlinear and include many different physical effects, it is not surprising to 
have multiple solutions. 

Consider as an example Fig. 2(b), which corresponds to an accretion luminosity log(Z acc ) = 
—0.75. The sequence of equilibria shows two peaks, one at log(/ ou t) ~ _ 2 and one at ~ —0.5. 
In addition, at the very right, there is a steep vertical segment which we refer to as the "wall." 
Fig. 4 shows other details of the equilibria for this choice of / acc ; panel (a) is a plot of Xi aycr 
vs the temperature T| ayer at the base of the layer, while panel (b) shows the variation of the 
H-fraction Xi aycr , He-fraction Fi aycr and heavy element fraction Z\ ayer at the bottom of the 
layer. 

When the column depth of the accreted layer £i ayer is very small, the gas is neither 
hot nor dense and does not undergo any nuclear burning. Thus, Xi ayer = X out , yi ayer = 
^out) -Ziayer = Z out , and the flux F out that emerges is mostly what is released as a result of 
compressing the accreting gas, plus any flux that escapes from the 10 8 K core. This stage 
corresponds to the nearly vertical segment at the left of Fig. 2(b). As £ increases, the gas 
at the bottom of the layer becomes hotter. Ultimately, when £ ~ 10 73 gem -2 and the 
temperature is about 10 7 ' 7 K, hydrogen-burning begins. Soon after this, hydrogen-burning 
takes over as the dominant term in the energy equation, and at this point, the £i ayer --F out 
curve in Fig. 2(b), as well as S layer -Ti ayer curve in Fig. 4(a), reverse direction, producing the 
peak on the left in the two plots. We identify this peak as the "hydrogen peak." The onset 
of hydrogen-burning is also evident in the run of Xi ayer in Fig. 4(b) (see also Fi ayer whose 
variations simply reflect the amount of helium produced by hydrogen-burning). 

As we descend from the top of the hydrogen peak in the direction of increasing F out and 
flayer, a point is reached when hydrogen-burning becomes saturated (because the reactions 
are beta-limited). This happens around T| ayer ~ 10 79 K. Beyond this point, the sequence of 
equilibria start rising again in the Si ayer -F out and £i ayer -7iayer planes. The rise continues for 
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Fig. 4. — Additional details of the equilibrium sequence shown in Fig. 2b. (a) Shows the 
column density S layer of the accreted layer as a function of the temperature T layer at the 
bottom of the layer. The open circles, the filled circles and the stars have the same meaning 
as in Fig. 2. (b) Shows the mass fractions of hydrogen, helium and metals at the bottom of 
the layer as a function of the normalized outgoing flux. 
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a while, with more and more hydrogen being burnt into helium until at Si ayer ~ 10 81 gem -2 
and Ti ay e r ~ 10 8 ' 3 K, helium-burning is initiated. At this point there is a second peak in the 
sequence of equilibria, the "helium peak." 

As the sequence of equilibria fall off from the helium peak, the helium is burned rapidly, 
and the CNO elements that this produces cause the hydrogen-burning also to pick up since 
the rate of hydrogen burning is proportional to Z CNO (eq. 26). When the hydrogen and 
helium are almost all exhausted, the sequence of equilibria turn round again and march up 
the wall. In this segment of the curve, very little changes as a function of increasing E. The 
hydrogen and helium are burned close to the surface within a column of order 10 8 gem -2 , 
and the rest of the accreted layer consists just of burned CNO material which is progressively 
compressed by the weight of the overlying gas. Note that, since helium burning by the triple- 
a reaction is not beta-limited, there is no regime of saturated helium burning. Therefore, the 
transition from the declining slope of the helium peak to the rapidly rising wall is generally 
quite abrupt. Figure 4(a) shows that, as Si ayer increases on the wall, the temperature at the 
base of the layer falls. This is because there is a net flux flowing from the layer into the star. 
The nuclear burning occurs at a fixed E in this sequence of wall models. With increasing 
Si ay e r , there is a larger dead column between the burning layer and the base of the layer, 
and the ingoing flux causes the temperature at the base to drop. 

We should note at this point that not all equlibria shown in Figs. 2-4 are accessible to a 
real system. Imagine starting with a bare neutron star and adding gas at the specified rate E. 
As Ei ayer increases, the system will ride up the left slope of the hydrogen peak. When Ei a y er is 
equal to the maximum E ~ 10 7 ' 3 g cm -2 of the hydrogen peak, the system no longer has any 
equilibria available in the vicinity of the peak. Therefore, it will relax towards the nearest 
available equilibrium, which is a solution with the same value of Si ay er ~ 10 7 ' 3 gem -2 on the 
left slope of the helium peak. Thus, all the equilibria on the right slope of the hydrogen peak 
and the lower left slope of the helium peak will be by-passed, since they correspond to lower 
values of Ei ay er than the current layer thickness. With increasing Ei a y er , the system will ride 
up the helium peak until it reaches the maximum E ~ 10 81 gem -2 of the helium peak. At 
this point, the system will once again move across to the nearest available solution, which 
is located on the wall, bypassing the equilibria in the valley to the right of the helium peak. 
Having reached the wall, the system will continue rising up the wall with increasing Ei a y er . 

Let us turn now to the other sequences of equilibria shown in Figs. 2 and 3. Fig. 2(a) 
corresponds to log(/ acc ) = —0.5. In this case, because of the rapid accretion rate, the accreted 
layer is quite hot, and so hydrogen burning is already in the saturated limit when it turns on. 
There is, therefore, no hydrogen peak (though there is the semblance of a plateau). Fig. 2(b) 
was discussed above, and has a modest hydrogen peak. Note that the helium peak in Fig. 
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2(b) is pushed to the right relative to the peak in Fig. 2(a). This is because helium requires 
a temperature of at least about 2 x 10 s K to burn (Fig. 1), and such high temperatures are 
obtained at the lower accretion rate only at larger values of F ont . With decreasing accretion 
rate, the gas temperature continues to fall and the hydrogen peak becomes more pronounced 
(Figs. 2c, d). The helium peak also gets pushed progressively farther to the right, until it 
finally hits the wall and disappears. For accretion luminosities below about 10~ L25 LEdd, 
there is no helium peak. Helium now starts burning on the wall, at progressively larger 
values of £i ayer with decreasing Z acc . 

The above results are for a core temperature of T core = 10 8 K. The pattern of peaks 
can be different for other core temperatures. If the core temperature is a few times 10 8 K, 
the hot core keeps the accreted layer hot under all circumstances and hydrogen burning is 
always in the saturated regime. In this case, there is no separate hydrogen peak. On the 
other hand, for lower values of T core , e.g. a few times 10 7 K, the hydrogen peak becomes 
quite pronounced and also has a higher peak value of Si ayer . In this case, it is often possible 
to have situations in which the helium peak is lower than the hydrogen peak. If such is the 
case, when the system reaches the top of the hydrogen peak it will move directly across to 
the wall without stopping at the helium peak. 



4. Stability 

The previous section described sequences of equilibria for various mass accretion rates. 
As explained in §2, for every equilibrium we carry out a linear stability analysis to check 
whether the system is stable or unstable (the latter is determined by the criterion given in eq. 
21). In Figs. 2 and 3, open circles indicate stable equilibria, filled circles indicate unstable 
equilibria in which the most unstable mode has a real 7 > 37 acc , and stars indicate unstable 
equilibria in which the fastest-growing mode has a complex 7, with ^(7) > 37 acc . 

From the various panels in Figs. 2 and 3, we quickly discern the following patterns. 
The left slopes of the hydrogen and helium peaks always correspond to stable equilibria, 
while the right slopes of both are always unstable; this is similar to the pattern described by 
Paczynski (1983a, who only considered He-burning) and is easily understood from his one- 
zone analysis. As explained in Appendix A, the same pattern is expected in our problem as 
well. The stability of equilibria on the wall is variable. Sometimes the equilibria on the wall 
are fully stable, sometimes they are fully unstable, and sometimes some equilibria are stable 
and some are unstable. As we discuss now, it is the stability of the wall that determines 
whether or not a particular accretion rate leads to thermonuclear bursts. 
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Consider as an example the case shown in Fig. 2(b), corresponding to log(/ acc ) = — 1. 
Starting from the bare neutron star, as gas piles up, the system moves up the hydrogen peak 
along a sequence of equilibria that are all stable (open circles). When Si ayer hits the top of 
the hydrogen peak, the system moves across to a point on the left slope of the helium peak 
with the same value of Si ayer - There is some adjustment of the accreted layer in order to 
switch to the new solution, but because the equilibrium solution here is stable the system is 
able to make the adjustment. With increasing £i ayer , the system climbs up the helium peak 
along a sequence of stable equilibria. As before, once it reaches the top of the peak, it moves 
across to the right, this time to the wall, and the solution here is again stable. However, as 
the system climbs up the wall, it reaches a point at which the equilibrium is unstable. At 
this point, there is no stable equilibria available for any value of /out- The system therefore 
undergoes a thermonuclear burst. 

We may now quickly identify which of the other cases shown in Fig. 2 and 3 are stable 
and which unstable. The high accretion rate system in Fig. 2(a) is stable because the states 
on the wall are all stable. The system will move across to the wall after it hits the helium 
peak, and it will then climb up the wall along a sequence of stable equilibria to arbitrarily 
large values of £i ayer - The cases shown in Figs. 2(b), 3(a),(b) are all similar. These systems 
reach the wall stably, climb part way up the wall and then go unstable; we call all these cases 
"delayed bursts." In contrast, the cases shown in Figs. 2(c), (d), 3(c), (d) become unstable 
the moment they hit the top of the last peak since there are no stable equilibria available on 
the wall; we call these "prompt bursts." 



4.1. Unstable vs Overstable Modes 

We pause to discuss the distinction between modes with real 7 and those with complex 7. 
The former correspond to a simple instability, where the mode amplitude grows exponentially 
with time, whereas the latter correspond to an overstability, where the mode oscillates at a 
frequency equal to the imaginary part of 7 even as the amplitude grows. As Figs. 2 and 3 
show, we see both kinds of behavior for our unstable modes. We should note that complex 
eigenvalues 7 always appear in complex conjugate pairs, which is obvious from the structure 
of the governing equations. 

As explained in Appendix A, from an inspection of the steady state and linear pertur- 
bation equations, it can be seen that any system that is at the top of a peak or the bottom of 
a valley, i.e., where c£i aycr / 'dF out = 0, has a zero-frequency mode, 7 = 0. Once we recognize 
this important rule, we can understand why the right slopes of the hydrogen and helium 
peaks are unstable. The argument goes as follows: the left slopes are stable, with all modes 



-22 - 



having < 0; the peak has a mode with 7 = 0; by continuity, this particular mode should 
have 7 > on the right slope and should be real. 

The above rule has a corollary: if the system makes a transition from stability to 
instability at a point where <i£i ayer /<iF out 7^ 0, then the unstable mode at this point must 
have a complex 7. This rule applies, for instance, to all the systems that have delayed 
bursts, i.e., become unstable half-way up the wall, e.g., Figs. 2(b), 3(a), (b). In these cases, 
when the system moves into the unstable zone, it always first hits an overstable region. 
The overstability may in some cases then become a pure instability as a pair of modes with 
complex 7 merge to spin off a pair of real 7. Also, it is not possible to tell "how complex" 
the mode is, i.e., the relative magnitudes of ^(7) and $5(7). 

Finally, in cases where the wall is entirely unstable, we have not been able to find any 
helpful rule. Whether the instability corresponds to real or complex 7 can be determined 
only with a full calculation. 

4.2. Burst Energetics and Trigger Mechanism 

Because there are two burning fuels, hydrogen and helium, we would like to know which 
fuel determines the burst properties. There are two aspects to this question. 

First, once an instability has been triggered, we would like to know how much energy 
is available for the burst in unburnt helium compared to the energy in unburnt hydrogen. 
Figure 5 shows this quantity as a function of accretion rate for three choices of the neutron 
star radius (assuming a mass of 1.4M Q ) and four core temperatures. We have assumed that 
all the energy in the unburnt fuel is emitted as radiation during the burst. In actuality, some 
energy is lost into the core and some comes out as neutrinos; we neglect this complication. 

In interpreting the results, we should note that hydrogen-burning to iron provides nearly 
5 times more energy per gram than helium-burning to iron. Keeping this in mind, we see that 
for high mass accretion rates roughly in the range log(/ acc ) ~ —0.6 to —1.5 (the precise values 
depend on the neutron star radius and core temperature), the bursts are of a mixed kind with 
energy being contributed by both hydrogen and helium. For intermediate accretion rates, 
log(^acc) ~ —1-5 to —2.5, we have pure helium-burning bursts. Finally, for log(/ acc ) < —2.5, 
we again have mixed bursts. Actually, this last range corresponds to pure hydrogen bursts 
and it is not clear if helium burning will be triggered during the burst. In making our 
estimates we assumed that all the unburnt fuel is consumed. The ranges of log(/ acc ) given 
here correspond to a core temperature of 10 s K and a neutron star radius of 10.4 km. 
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Fig. 5. — Ratio of the fluence in the burst due to helium burning to the fluence from 
hydrogen burning, assuming that each fuel is burned to 56 Fe. The four panels correspond 
to core temperatures T core = 10 8 ' 5 , 10 8 , 10 7 ' 5 and 10 7 K, respectively. In each panel, the 
triangles, circles and stars trace the results for a neutron star of mass 1.4M Q and radii 
R = 10 (i Rs = 16.4 km, lO - 4 ^ = 10.4 km, and lO - 2 ^ = 6.5 km, respectively. Note that 
hydrogen and helium burning both contribute at high and low luminosities (with hydrogen 
being about five times more important energetically per unit mass), while helium burning 
dominates at intermediate luminosities log(/ acc ) ~ —1.5 to —2.5. 
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Another interesting question is to understand what precisely triggers the burst insta- 
bility. One way to do this is to look at the maximum temperature in the fuel layer (Fig. 6) 
and to thereby identify which fuel might be important. As we discussed in §2.3, hydrogen 
burning becomes saturated, and hence stable, for temperatures above about 10 7 9 K (see 
Fig. 1). Therefore, we expect hydrogen-triggered bursts at lower temperatures and helium- 
triggered bursts at higher temperatures. From Fig. 6 we see that the former occur at lower 
accretion rates and the latter at higher accretion rates. However, there is no clear way to 
distinguish the kinds of bursts from this plot since the variation of T max with accretion lu- 
minosity is smooth. In our models, we do not find unstable modes for layers with maximum 
temperatures above 3.5 x 10 8 K. This temperature limit is reached at log(/ acc ) ~ —0.6. 

Gaining a deeper understanding of what triggers the bursts requires an examination 
of the eigenf unctions of the unstable modes. Specifically, we compute the contribution of 
helium burning integrated over the eigenfunction and take the ratio of this to the hydrogen 
burning integrated over the eigenfunction (here we do not consider the additional burning 
to iron as we did with the burst fluence ratio). The ratio should be large for a pure helium- 
triggered burst and small for a pure hydrogen-triggered burst, and on the order of 0.1 — 0.2 
for mixed-triggered burst (the ratio is not 0.5 in this case since hydrogen-burning releases 
nearly an order of magnitude more energy per gram than helium-burning). 

Figure 7 shows the above ratio as a function of log(/ acc ) for the same three choices of 
the neutron star radius and four choices of the core temperature. We see that, at high 
mass accretion rates log(/ acc ) ~ —0.6 to —1.5, the bursts are triggered by both hydrogen 
and helium. Within this range, the higher accretion rates give delayed bursts, while the 
lower rates give prompt bursts (as seen from Fig. 2). For intermediate accretion rates 
log(/ acc ) ~ —1-5 to —2.5, we have pure helium-triggered bursts. These all involve overstable 
modes. Finally, for log(/ acc ) < —2.5, we have hydrogen-triggered bursts. These are pure 
unstable modes (real 7). 

Figures 8-9 show a typical eigenmode corresponding to log(/ acc ) = —0.9. Near the sur- 
face (£ — > 0), all perturbations vanish according to our boundary conditions, except the flux, 
which is given a unit perturbation (in arbitrary units — it merely sets the normalization). 
The accreted layer extends up to Si ayer , which is indicated by the thick vertical line. The 
eigenfunction then continues on into the substrate to a depth such that the diffusion time to 
that point is equal to twice the mode time scale (see §2.2). Note that, even though the inner 
boundary condition corresponds to a vanishing temperature perturbation, in fact the flux 
perturbation also vanishes at this point. This shows that our choice of the inner boundary 
is physically well-motivated. We have tried integrating the eigenmode to larger or smaller 
depth in the substrate and generally obtained the same results. 
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Fig. 6. — The maximum temperature T max achieved in the accreted layer when the layer first 
becomes unstable, plotted as a function of the accretion luminosity. The core temperature 
is assumed to be T core = 10 8 K, and the calculations correspond to the equilibrium at the 
"wall." The symbols have the same meanings as in Fig. 5. 
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Fig. 7. — Ratio of the energy produced by helium burning in the linear eigenmode to the 
energy produced by hydrogen burning. The symbols and the layout are the same as in 
Fig. 5. Note that at high luminosities there is mixed burning with significant contributions 
from both hydrogen and helium, at intermediate luminosities helium burning dominates, 
while at low luminosities hydrogen burning dominates. 
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Fig. 8. — Overstable eigenmode for log(/ acc ) = —0.9, M = 1.4M , R = 10.4 km and T core = 
10 8 K. The solid and dashed lines trace the real and imaginary parts of the mode. The bold 
vertical line shows the bottom of the accreted layer at £ = £i aycr . The eigenmode continues 
into the substrate to a depth £ max , at which point the amplitude of the mode is substantially 
decreased. The boundary condition at £ max is that the temperature perturbation should 
vanish. However, it is seen that the perturbations in essentially all variables vanish there. 
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Fig. 9. — Expanded version of the eigenmode in Fig. 8, showing the burning layer in more 
detail. 
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Within the accreted layer, the finite flux perturbation at the surface causes a tempera- 
ture perturbation over the entire layer. However, the perturbations in X and Z are limited 
to a very narrow layer. This is the burning layer, which is shown in an expanded view in 
Fig. 9. The fluctuation in X indicates hydrogen burning and the fluctuation in Z corresponds 
to helium burning. For this mode, we see that both fuels burn at roughly the same depth, 
indicating that it is an example of a mixed burst. Also, since the eigenmode has nontrivial 
real and imaginary parts, we have a complex eigenvalue 7. In this particular case, the real 
part of 7 is equal to 8.9 x 1CT 5 s _1 and the imaginary part is equal to 6.6 x 10~ 4 s _1 . 

At lower accretion rates, the eigenmodes tend to be simpler. If we consider the critical 
values of Si ayer at which the systems shown in Figs. 3(a),(b) become unstable, the burning 
is dominated by helium and there is only a small contribution from hydrogen. Moreover, the 
two fuels tend to burn at different values of S. The eigenmode is still complex, though the 
imaginary part is not very large. If we go to yet lower accretion rates as in Figs. 3(c), (d), 
then there is no helium burning at all and the burning layer is dominated by hydrogen. In 
this case, the eigenmode is fully real. 



4.3. Putting Things Together 

Pulling together the various ideas discussed above, we identify five regimes of bursting 
behavior as a function of decreasing accretion luminosity. In the following, the numerical 
values correspond to the particular calculations presented earlier, in which M = 1.4M , 
R = 10 0A R S = 10.4 km, T core = 10 8 K. Tables 2 and 3 give results for other selected choices 
of R and T core . 

I. For very high accretion rates, with log(/ acc ) > —0.6, there are no thermonuclear bursts. 
This is the regime of stable accretion, where the accreting gas is able to burn the nuclear 
fuel without instability. 

II. For somewhat lower accretion rates log(7 acc ) ~ —0.6 to —0.85, we have mixed bursts 
triggered by both hydrogen and helium. These bursts burn substantial amounts of both 
fuels. In this particular accretion range, the instability is delayed and is triggered only after 
the system climbs part way up the wall. We call these "delayed mixed bursts." As the system 
climbs up the wall, a large fraction of the nuclear fuel is burned stably in the accreted layer 
before the burst is triggered. 4 As we shall see, this has dramatic consequences for both the 



4 Note that, on the wall, / ou t is nearly equal to unity, which means that the escaping flux in the equi- 
librium solution is almost equal to the nuclear energy content of the new fuel being added to the surface. 
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recurrence times of the bursts and the ratio of the burst fluence to the total accretion fluence 
between bursts. When the system becomes unstable, the mode has a complex 7, which 
corresponds to an overstability. Thus, we expect the system to show an oscillatory behavior 
as it approaches the instability. 

III. For yet lower accretion rates log(/ acc ) ~ —0.9 to —1.25, we continue to have mixed 
triggered bursts burning mixed fuel. However, here the entire wall is unstable and so the 
instability is triggered when the system reaches the top of the helium peak. Because the 
burst happens as soon as the system hits the wall, we refer to these as "prompt mixed 
bursts," to distinguish them from the previous class of delayed mixed bursts. When the 
instability first begins at the peak, 7 is real, but when the system reaches the wall, 7 may 
possibly be complex (though not always). We cannot tell from our present analysis whether 
or not the system will show oscillatory behavior prior to the burst. 

IV. For still lower accretion rates log(/ acc ) ~ —1.3 to -2.5, bursts are again triggered only after 
the system climbs part way up the wall, i.e., we have delayed bursts. Here the instability 
is helium-triggered and the burst is dominated by helium burning. These correspond to 
the classic "helium bursts" that have been studied by previous authors. According to our 
calculations, 7 is complex when these systems go unstable, and this means that there should 
be some kind of oscillatory behavior prior to the burst. In practice, the imaginary parts of 
7 tend to be small for these modes. 

V. Finally, for the lowest accretion rates log(/ acc ) < —2.5, the wall is again unstable, and now 
the burst happens as soon as the system hits the top of the hydrogen peak. The instability 
is triggered by hydrogen burning, and we call these "hydrogen bursts." They do not have 
any oscillatory behavior. It is not clear if these hydrogen bursts will trigger helium burning. 
If they do not, then the helium will accumulate as the ashes of hydrogen burning and the 
system is likely to have enormous but very rare helium bursts. 

The various regimes described above are summarized in Tables 2, 3. Much of our 
discussion in the rest of the paper is in terms of the burst regimes discussed in this section. 
Many of the regimes have been recognized by previous authors (Fujimoto et al. 1981; Fushiki 
& Lamb 1987a; Bildsten 2000), but one of the new ideas to come out of our analysis is the 
distinction between delayed mixed bursts and prompt mixed bursts. This distinction has 
not been made earlier. The other new aspect is our ability to distinguish between simple 
exponential instability and oscillatory overstability. It is unique to our work since this is the 



Nevertheless, the accreted layer does possess energy in unburnt fuel. Loosely speaking, it is the unburnt 
fuel that accumulated in the layer prior to reaching the wall. The energy from burning this fueld is released 
during the burst. 
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Table 2. Burst Regimes for T core = 10 7 5 K 







Range in 


log(/ acc ) 










R = 6.5 km 


R = 10.4 km 


Type of Burst 


Symbol in Fig. 10 


I. 


> 


(-0.7) 


> (-0.55) 


No Bursts 


None 


II. 


(- 


-0.75) to (-1.0) 


(-0.6) to (-0.85) 


Delayed Mixed 


Naked Star 


III. 


(- 


-1.05) to (-1.6) 


(-0.9) to (-1.4) 


Prompt Mixed 


Circled Symbol 


IV. 


(- 


-1.65) to (-2.3) 


(-1.45) to (-2.4) 


(Delayed) Helium 


Naked Star 


V. 


< 


(-2.35) 


< (-2.45) 


(Prompt) Hydrogen 


Circled Symbol 



Table 3. Burst Regimes for T corc = 10 8 K 







Range in 


log(/ acc ) 










R = 6.5 km 


R = 10.4 km 


Type of Burst 


Symbol in Fig. 10 


I. 


> 


(-0.7) 


> (-0.55) 


No Bursts 


None 


II. 


(- 


-0.75) to (-1.0) 


(-0.6) to (-0.85) 


Delayed Mixed 


Naked Star 


III. 


(- 


-1.05) to (-1.4) 


(-0.9) to (-1.25) 


Prompt Mixed 


Circled Symbol 


IV. 


(- 


-1.45) to (-2.4) 


(-1.3) to (-2.5) 


(Delayed) Helium 


Naked Star 


V. 


< 


(-2.65) 


< (-2.55) 


(Prompt) Hydrogen 


Circled Symbol 
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first study to carry out a full linear stability analysis to calculate eigenmodes and (complex) 
eigenfrequencies. 

The other comment we should make is that the pattern of five regimes described above 
is valid only for intermediate core temperatures ~ 10 7 5 — 10 8 K. We have focused on this 
case since it represents the likely core temperature for standard neutrino cooling models of 
accreting neutron stars. For higher core temperatures, e.g., 10 85 K, the hydrogen-burning 
layer is hotter than 10 78 K so that the systems are in the saturated hydrogen-burning limit. 
They therefore lack the hydrogen-triggered bursts at low accretion rates. On the other hand, 
if the core temperature is lower, e.g., 10 7 K, the prompt mixed bursts are no longer triggered 
by the helium peak. Instead, the bursts are triggered when the accreted layer reaches the 
top of the hydrogen peak, which is now higher than the helium peak. The bursts are still 
of the mixed variety though. In addition, there are quantitative differences as a function of 
neutron star radius and core temperature. These are discussed in the next section. 



5. Results 

5.1. Different Regimes of Bursting 

Figure 10 depicts the various bursting regimes as a function of the radius of the neutron 
star and the accretion rate. The results shown are for core temperatures of T core = 10 8 K and 
10 7 K. The blank areas correspond to stable nuclear burning, the symbols enclosed within 
circles correspond to prompt bursts and the naked symbols correspond to delayed bursts. 
The five regimes of bursting behavior described in §4.3 and summarized in Tables 2, 3 are 
clearly delineated. 

Moving from right to left in the panel corresponding to 10 s K, we see that at the highest 
mass accretion rates there is no burst activity since these systems consume hydrogen and 
helium stably (region I in §4.3). Immediately to the left is a region of delayed mixed bursts 
(region II). At lower accretion rates, there is a zone of prompt mixed bursts (symbols with 
circles, region III). Then there is a relatively broad zone of delayed helium bursts (region 
IV), and finally, at the lowest accretion rates, a region of prompt hydrogen bursts (symbols 
with circles, region V). 

Some of the boundaries between the different zones are remarkably insensitive to the 
neutron star radius, though the zone of prompt mixed bursts does vary with radius, at 
least for T core = 10 s K. For very compact neutron stars, e.g., \og(R/ R s ) = 0.2, this zone 
is reasonably wide, whereas for \og(R/R s ) = 0.6, the zone almost disappears. The zone of 
prompt mixed bursts is much broader for T core = 10 7 K. These variations with radius and 
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T core lead to noticeable signatures in many of the diagnostics that we discuss later. 

In Fig. 10, the points represented with stars correspond to overstable modes with com- 
plex 7 and the triangles correspond to unstable modes with real 7. In accordance with the 
discussion in §4.3, we see that delayed bursts are nearly always overstable, while prompt 
bursts are either overstable or unstable. 

The map of the different burst regimes is roughly the same at other core temperatures, 
though the widths and locations of the various zones tend to vary. Generally, the regime of 
delayed mixed bursts becomes somewhat narrower with decreasing T core . Also, for T corc = 
10 8 ' 5 K, there are no hydrogen bursts, i.e., the leftmost region of filled triangles within circles 
is absent. These patterns may be noticed in several of the plots discussed later. 



5.2. Burst Recurrence Time 

We begin our overview of the observable properties of thermonuclear bursts on neutron 
stars with the recurrence time t rec , i.e., the average time between bursts. Because we only 
treat the physics of the accreted layer until a burst is initiated, but not the burst itself, we 
need to make some assumption as to how much of the available fuel is consumed in a burst. 
In our calculations we assume that all the fuel is burned, which means that our estimates 
of t rec are really upper limits to the observed recurrence time. With this approximation, if 
Siaycr.crit is the critical column depth of the accreted layer at which a burst is initiated, then 
the recurrence time t rec measured by a distant observer is given by 

t rcc =(l + z)^^. (33) 

Figure 11 shows t rcc as a function of the accretion luminosity for three choices of the 
neutron star radius, \og(R/ Rg) = 0.2, 0.4, 0.6 (i.e., R = 16 A, 10.4, 6.5 km), and four core 
temperatures, T core = 10 8 ' 5 , 10 8 , 10 7 ' 5 , 10 7 K. The 10 8 K panel, and to some extent the 10 7,5 
K panel, clearly show the different regimes of bursting behavior discussed in the previous 
subsection, especially for smaller neutron star radii. The recurrence time tends to become 
larger for delayed bursts compared to prompt bursts. The straight band of points between 
log(Zacc) ~ —1 and —1.6 for the two smaller radii corresponds to prompt mixed bursts. This 
band is missing for the largest radius because, as seen in Fig. 10, the region of prompt mixed 
bursts disappears for that radius. 

One obvious pattern in Fig. 11 is that there is a general overall increase of t rec with 
decreasing L acc . This is because £ becomes lower with decreasing accretion rate, so it takes 
longer to build up a given accretion column. To illustrate this point, Fig. 12 shows the 
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Fig. 10. — Regimes of burning on the surface of a 1.4M neutron star for T core = 10 8 K and 
10 7 K. The empty areas correspond to stable burning, the filled triangles denote regions with 
unstable modes (real 7), and the stars to regions with overstable/oscillatory modes (complex 
7). Prompt bursts (see the text for an explanation) are shown by symbols circumscribed 
with a circle, and delayed bursts are shown by symbols without circles. Starting from the 
right, with decreasing luminosity, one goes through the following sequence: stable, delayed 
mixed bursts, prompt mixed bursts, helium bursts, hydrogen bursts (compare with Table 3). 
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Fig. 11. — The burst recurrence time t rec for a distant observer, assuming that all the fuel is 
consumed. The four panels correspond to four core temperatures, and the triangles, circles 
and stars correspond to neutron star radii of 16.4 km, 10.4 km and 6.5 km, respectively. The 
four regimes of bursting described in §4.3 and summarized in Tables 2, 3 are clearly seen in 
many of the sequences. 
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Fig. 12. — Critical column density of accreted material S layer crit at the onset of the burst. 
The symbols and layout are the same as in Fig. 11. 
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dependence of £i ay er,crit for the various cases. We see that, for the prompt mixed bursts, 
Diayer.crit is more or less constant. 

The results for the other core temperatures are generally reasonable. For T corc = 10 8 5 K, 
the recurrence times are fairly similar to the 10 8 K case, except that there are no hydrogen 
bursts at low mass accretion rates. The case of T core = 10 7 - 5 K is very similar to 10 s K. 
However, for T corc = 10 7 K, the column needed to produce bursts tends to be higher. This 
leads to significantly longer recurrence times. 

The shortest recurrence time we find in the model is about 5 hours, which occurs at 
log(Zacc) ~ — 1 for intermediate core temperatures T core ~ 10 7 5 — 10 s K. The corresponding 
maximum burst rate is about 5 bursts/day. 



Probably the most easily observed property of bursts is the total fluence, Eh + E He , 
resulting from burning the fuel in the accreted layer. We assume that all of the unburnt 
hydrogen and helium in the accreted layer is fully burned to iron in the burst and radiated 
during the burst (neglecting any flux going into the star or any energy loss through neutrino 
emission). We then divide the burst fluence by the Eddington luminosity to calculate an 
effective burst duration: 



This quantity is plotted in Fig. 13 for our standard three neutron star radii and four core 
temperatures. As in the case of the recurrence times discussed above, the calculated values 
are upper limits because of the assumption of complete burning of the fuel. 

Since hydrogen burning is typically beta-limited during the burst, the hydrogen fluence 
may in some cases emerge as a long plateau following the initial spike of helium burning 
(which is not beta-limited). Therefore, another possibly useful measure of burst duration is 
the effective duration t^e of helium burning alone, calculated again as fluence divided by the 
Eddington luminosity. This quantity is plotted in Fig. 14. 

Consider first the regime of helium bursts, which corresponds to the segments of tH+He 
and £hc with steep negative slopes in the range log(/ acc ) ~ —1.5 to —2.5. These are delayed 
bursts in which the hydrogen burns stably at a shallow depth, while the helium accumulates 
over a fairly large column before the system goes unstable. The bursts are very much 
dominated by helium burning (see Fig. 5), and so tH+He and tu e are nearly equal. Also, the 
fluence is directly proportional to the accretion column Si aycr , as can be seen by comparison 
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Fig. 13. — Effective burst duration, calculated by dividing the total burst ffuence by the 
Eddington luminosity. Triangles, circles and stars correspond to neutron star radii of 16.4 km, 
10.4 km and 6.5 km respectively. Notice that, at the extreme right of the various panels, 
where we have the regime of delayed mixed bursts, the burst duration is either level or falls 
relative to the prompt mixed bursts, even though the recurrence time t rec (Fig. 11) and 
the layer thickness Si ayer crit (Fig. 12) increase substantially. This indicates that there is 
considerable steady burning of nuclear fuel in this regime. 
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with Fig. 12. 

There are two segments of prompt bursts, one at higher luminosities and the other at 
lower luminosities, on either side of the helium bursts. Both involve significant unburnt 
hydrogen and helium, and therefore for both tH+He is significantly greater than tu c (since 
hydrogen burning releases more energy per gram, by a factor of about 5, compared to 
helium burning). Both burst durations are roughly proportional to Xi aycr as can be seen by 
comparison with Fig. 12. 

Finally, the delayed mixed bursts, which are present for log(7 acc ) > — 1, show a different 
behavior. These are systems in which considerable stable burning of nuclear fuel occurs in 
the accreted layer before the layer actually bursts. Therefore, the amount of nuclear fuel 
available to power the burst is much less than one might expect for the given column depth. 
The burst durations are at most comparable to, and are often less than, the durations of 
adjacent prompt bursts, even though these delayed bursts have larger values of £i aycr and 
correspondingly much longer recurrence times (compare Figs. 13, 14 with Fig. 12). 

From Figs. 13, 14 we see that the typical burst durations are predicted to be a few 
seconds to few tens of seconds at high mass accretion rates. However, the duration can be 
quite large, an hour or longer, for helium bursts with log(/ acc ) < —2, and as much as a day 
for cooler core temperatures. 



5.4. Burst a 

The recurrence times and burst durations discussed in the previous subsections were 
calculated assuming that the entire fuel reservoir is consumed in a burst. However, this 
assumption may not be valid (e.g. Taam 1982; Fujimoto et al. 1987b), especially for cold 
neutron stars (Taam et al. 1993). It is, therefore, useful to consider the ratio a of the nuclear- 
burning energy that is emitted during the burst to the total energy released between bursts 
by accretion. In terms of quantities we have introduced earlier, a is given by 

t rec Z/ acc t rec /or\ 

« = f _i_ f = 7 Z aoc- (35) 

-^H + -^He f-H+He 

The quantity a should be independent of what fraction of the fuel is burned in bursts; if only 
a fraction of the fuel is burned, then both t rec and £h+Hc will decrease by the same fraction 
and a ought to remain the same. This makes it a particularly useful parameter. 

Fig. 15 shows a from our models for several neutron star radii and core temperatures, 
and a range of accretion rates. Two trends are obvious. First, more compact neutron stars 
have larger values of a compared to less compact stars; the former release more gravitational 
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Fig. 14. — Similar to Fig. 13, but for only helium burning. 
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energy per unit accreted mass compared to the latter, whereas both release roughly the same 
amount of nuclear energy in bursts. Second, prompt bursts tend to have smaller values of 
a compared to delayed bursts. Generally, prompt bursts have a lot of unburnt hydrogen 
available. This gives a larger burst fluence per unit mass and thus a smaller a. The delayed 
helium bursts have very little hydrogen. Since helium burning releases a factor of about 5 
less energy per gram, the a of these bursts is larger by about this factor. The increase of a 
for the delayed mixed bursts at log(/ acc ) > —1 is for a different reason. Here, the accreting 
layer has both hydrogen and helium, but both fuels are burned substantially before the burst 
is initiated. This is the reason for the upturn of a at the highest accretion rates in all the 
panels in Fig. 15. This trend has not been noted in previous work. 



5.5. Constraints on Neutrino Cooling 

If we assume that the neutron star is in thermal equilibrium as it accretes we may 
estimate the temperature of its core as a function of its mean accretion rate. There are 
two contributions to energy flowing into the core. The first is the inward flux at the inner 
boundary of the accreted layer and the substrate; this may be either positive or negative 
(see Brown 2000). The second is the energy of released by nuclear reactions deep in the 
crust (Brown et al. 1998), typically ~ 1 MeV/baryon. Our calculations show that the 
former contribution is much smaller than the latter; this is consistent with the discussion in 
Cumming & Bildsten (2000) where the authors estimate that, out of the 1 MeV per baryon 
released by deep crustal nuclear reactions, only about 150 keV per baryon escapes to the 
surface. Therefore, for an approximate estimate of the core temperature, it is reasonable 
to assume that the flux flowing into the core is equal to 1 MeV/baryon (see Brown 2000 
for a more detailed discussion including the effect of different assumptions on the thermal 
conductivity). If we further assume the simple estimates of the neutrino emissivity given in 
Shapiro & Teukolsky (1983), we may estimate the equilibrium temperature of the core for a 
given accretion rate. For the modified URCA process we find 

^y_ 2500 (i+iii^-f-g.^, (36) 



y nuc 



10 s KJ ~ 1 MeV/baryon z M Edd \pn 

which gives a typical temperature in the range (1 — 3) x 10 8 K. If the neutron star cools via 
the direct URCA process from (say) a pion core, then the core temperature is 

r ~ ^ - 9 £ °" V + *?_*_ JLr«, (37) 



10 7 KJ IMeV/baryon z M Edd p nuc 

where 9 ~ 0.3. In this case, the core is significantly cooler, T corc ~ 10 7 K. In both cases, 
the core temperature depends only weakly on the mean accretion rate (as M 1//§ and M 1//6 , 
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Fig. 15. — Dependence of the parameter a, the ratio of the energy generated by accretion 
between bursts to the nuclear energy emitted during bursts. Triangles, circles and stars 
correspond to neutron star radii of 16.4 km, 10.4 km and 6.5 km, respectively. Notice the 
steep rise of a at high accretion luminosities, near the stable burning limit. The rise is 
characteristic of delayed mixed bursts. 
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respectively) and on the exact flux flowing into the core (whether it is 1 MeV per baryon or 
a little more or less.) 

An important caveat is that our results for cool neutron stars with T corc ~ 10 7 K depend 
somewhat sensitively on how we treat the inner boundary condition (see §2.2), so the bursting 
behavior of stars with direct URCA cooling may not be precisely as depicted by our results 
for 10 7 K; however, it should be noticeably different from that of stars with modified URCA 
cooling. 

6. Discussion 

We begin in §6.1 with a brief summary of how the results depend on the input physics. 
We follow that in §6.2 with a discussion of different kinds of modes present in our model 
and use this to explain the occurrence of delayed mixed bursts. In §6.3, we discuss how our 
methods compare with the many previous theoretical analyses of Type I bursts in neutron 
stars. We then present in §6.4 a preliminary comparison of our theoretical predictions with 
observations, making extensive use of the review of EXOSAT data published by van Paradijs 
et al. (1988). We defer a detailed comparision with more recent observations with the Rossi 
X-ray Timing Explorer and BeppoSAX to a subsequent paper, but briefly touch on the work 
of Cornelisse et al. (2003). (The latter paper was posted after our paper was submitted 
to the journal. Since their results are very relevant for this work, we have included some 
discussion.) Finally, in §6.5, we discuss some improvements to this work that may be worth 
pursuing. 

6.1. Sensitivity to Input Physics 

In the course of doing this work, we changed several prescriptions for the input physics, 
usually starting with a simple approximation and graduating later to more realistic prescrip- 
tions. By monitoring how the results changed we have developed a sense of which aspects 
of the input physics are most important if one is interested in accurate results. 

In the equation of state, the only complication is the electron pressure. Before settling 
on the quadrature prescription given in equation (22), which is taken from Paczynski (1983b), 
we tried a simpler prescription in which we wrote the electron pressure as the straight sum 
°f -Pc,nd and P e ,d- The differences in the results are insignificant. Similarly, we tried modeling 
the nuclear reactions both with and without screening. There was little difference, except 
perhaps at the lowest accretion rates, log(/ acc ) ~ —3, where the hydrogen bursts showed 
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some small variations. 

The effect of changing the opacity prescriptions is more serious. For instance, the 
opacities of Hernquist & Applegate (1984), which we tried first, are generally lower than 
those of Potekhin (1999). Therefore, with the former, a thicker layer must accumulate 
before the layer becomes unstable. This effect is strongest for the hydrogen bursts, and 
yields longer recurrence times and larger fluences. Also, because the burning layer is more 
poorly insulated from the core, the results tend to be more sensitive to the core temperature 
than when we use the modern conductive opacities of Potekhin (1999). As a result, with the 
Hernquist & Applegate (1984) opacities we find that models with T corc = 10 7 ' 5 K are more 
similar to models with T core = 10 7 K, whereas with the Potekhin (1999) opacities, they are 
more similar to 10 8 K. 

The results are even more sensitive to the radiative opacity prescription. Prior to 
settling on the Schatz et al. (1999) prescription, we employed the simple fitting function 
of Iben (1975). The results were qualitatively the same as the ones presented here, except 
that the boundaries between the different regimes moved towards lower values of log(/ acc ) by 
about 0.2. That is, the patterns in Figs. 5, 10-15 were shifted to the left by this amount, 
but not altered very much in shape. In view of this, it may be worthwhile to incorporate an 
even better approximation for the radiative opacity in future calcualtions. 

The frequencies and structure of the unstable modes are not changed significantly by 
the opacity with the exception that the modes penetrate deeper into the substrate for lower 
opacities. 

Finally, we have found that the inner boundary condition plays an important role. 
Figure 16(a) shows two sets of calculations, both of which set the temperature equal to T core 
at the inner boundary. In one case, the matching is done inside the substrate at a depth E di g 
determined by a diffusion criterion (see §2.2 and eq. 18 for details), while in the other the 
boundary condition is set at Si ayer , the bottom of the accreted layer. There is an enormous 
difference in the results. In the former case, which is more physical, the temperature of the 
accreted layer tends to remain high and the temperature profile relaxes to T corc only well 
inside the star. In the latter case, however, the boundary condition forces the gas layer to be 
cool near the bottom and this causes the rest of the layer to be cooler than in the previous 
case. The lower temperature delays the onset of nuclear burning, causing the critical surface 
density for instability to go up significantly. 

We feel that we have captured a good fraction of the important physics with the bound- 
ary condition that we use. Nevertheless, there is probably room for further improvement. 
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Fig. 16.— Results for T corc = 10 8 K, R = 10.4 km, M = 1.4M . (a) Shows how the 
restriction to thermal modes (filled circles) affects the results. Compared to the full instability 
calculation (open circles), only a fraction of the range of / acc is found to be unstable for purely 
thermal perturbations. The results depicted by crosses use the instability criterion of Fushiki 
& Lamb (1987a). The lower set of circles corresponds to our standard boundary condition 
in which the temperature boundary condition is applied inside the star at a diffusion depth. 
The upper set of circles depicts the results if the temperature boundary condition is applied 
at the bottom of the fuel layer [as in Narayan & Heyl (2002) and Fushiki & Lamb (1987a)]. 
(b) Shows how the choice of g mo d e affects the onset of instability. The larger the value of 
fl'mode, the more the nuclear modes near the right end of the plot are delayed. The thermal 
modes on the left are hardly affected. 
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6.2. Thermal Modes, Nuclear Modes, and Delayed Mixed Bursts 

Because the eigenmodes and growth rates are calculated by perturbing five coupled 
differential equations, the physics of the modes is not always apparent. We make a beginning 
here by sorting the modes into two major classes. 

To understand the division, we note that the perturbation equations written in Appendix 
A involve three time derivatives, one each in the energy equation (A10), the H-evolution 
equation (A6) and the He-evolution equation (A12). The characteristic time scale of the 
energy equation is something like the diffusion time, which is usually quite short, whereas 
the characteristic time scale of the other two equations is the accretion time, since this is the 
time on which the composition of a parcel of gas changes. We thus expect fast-growing modes 
to mostly involve the the energy equation (plus the equations of hydrostatic equilibrium and 
energy transfer which do not involve time derivatives), while slower modes should involve 
significant perturbations in all variables, including X and Y . We refer to the former as 
"thermal modes" and the latter as "nuclear modes." 

It is straightforward to filter out the nuclear modes in the calculations — we just switch 
off perturbations in X and Y by setting X(Y0) = X (S), Y(E) = Y (E), or equivalently, 
Xi(E) = Y]_(E) = 0. Figure 16(a) shows the results of such a calculation. Compared to our 
standard calculation in which all variables are allowed to vary (open circles), we see that 
thermal modes (filled circles) are present only for lower accretion rates, log(/ acc ) < —1.5, i.e., 
in the regime of helium bursts and hydrogen bursts, but not at higher accretion rates where 
mixed bursts occur. An inspection of the mode growth rates confirms that these modes grow 
on relatively short time scales; typically, the growth time for thermal modes is a factor of 
tens to hundreds shorter than the accretion time scale. The modes above log(Z acc ) ~ —1.5, 
which are present only when all variables are allowed to vary, are the nuclear modes in our 
classification. These modes have slow time scales, and their eigenfunctions involve important 
fluctuations in composition, which is not the case with the thermal modes. The eigenfunction 
shown in Figs. 8, 9 is a nuclear mode. 

The identification of the slow nuclear modes provides a natural explanation for the 
category of delayed mixed bursts that we have newly identified in this paper. Recall that we 
do not consider a system to be unstable unless the growth rate ^(7) is greater than g mo de7acc 
(see eq. 21). For our nominal choice of g mo d e — 3, the thermal modes are not in the least 
affected by the value of gmode, since their growth rates are typically much greater than the 
limit. However, the slower nuclear modes are strongly influenced by the criterion, and the 
effect is particularly severe at higher accretion rates. The right panel of Figure 16 shows the 
critical column for bursts Si ayer crit as a function of log(/ acc ) for different choices of <? mo de- The 
delayed burst regime becomes more and more prominent as g mo dc increases. The calculations 
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show that, for high accretion rates, the accreted layer is marginally unstable already at small 
values of S layer . However, the growth rate is low, and the instability takes a while to grow. 
By the time the mode has grown enough to produce a burst, S laycr is significantly larger; 
this is the reason for the delay in the burst. In addition, because of the delay, much of the 
fuel is burned stably prior to the burst, and the amount of unburnt fuel available for the 
burst is reduced. This causes a dramatic increase of a in this regime. The extent of these 
effects depends on the choice of g mo d e - The value we have used for the results presented in 
§5, .^mode = 3, are in our opinion reasonable. 

6.3. Comparison with Earlier Theoretical Work 

Previous investigations have identified the various burning regimes that we have found 
here (Fujimoto et al. 1981; Fushiki & Lamb 1987a; Bildsten 2000), with the exception of the 
regime of delayed mixed bursts (see §4.3, Tables 2, 3, Fig. 10). The ranges of Z acc of the 
different regimes obtained by these workers generally agree with our results, though we have 
found several new features, as discussed in previous sections. Our method, being more rig- 
orous, promises to provide better quantitative predictions for comparison with observations. 
Also, we calculate mode frequencies, which enables us to identify overstable modes and to 
estimate the oscillation periods. It is unclear if this will have clear observational signatures, 
but it is a topic worthy of further investigation. 

As described in §1, our method is similar in spirit to previous investigations that eval- 
uated steady-state configurations and then studied the stability of the equilibria to small 
perturbations (for example Hansen & van Horn 1975; Ergma & Tutukov 1980; Fujimoto 
et al. 1981; Fushiki & Lamb 1987a; Cumming & Bildsten 2000). The computation of the 
equilibria in the various studies probably do not differ a great deal since the physics is ba- 
sically well understood (§2). The approximations that we and others use (§2.3) appear to 
be harmless, except perhaps the opacity which does make a difference (§6.1). Once one 
has calculated a sequence of equilibria such as those shown in Figs. 2, 3, it is necessary to 
determine at what £i aye r,crit, if any, the accreted layer becomes unstable. It is at this stage 
that the major differences in methods appears. 

The technique that we have used in this paper is new and avoids any prejudices as to 
which processes are important and which are not for triggering a burst. We carry out a full 
linear stability analysis of the equilibrium solution by solving for the eigenmodes and their 
complex eigenfrequencies. In particular, the regime of delayed mixed bursts that we have 
found (§4.3, Tables 2, 3), and the identification of nuclear modes (§6.2), are entirely the result 
of the more rigorous formalism we employ. As we shall see below, the delayed mixed burst 
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regime may help to explain some puzzling observations. Also, none of the previous studies 
was able to identify whether the growing mode is a simple instability or an overstability. 
This too may have some observational implications. 

Earlier works have generally used various local or approximately global criteria, calcu- 
lated at or near the bottom of the layer, to determine the stability of the accreted layer. 
With the exception of Fushiki & Lamb (1987a), no one to our knowledge has discussed global 
perturbations of the steady-state configuration. Even the Fushiki & Lamb (1987a) study was 
relatively crude since it assumed a constant temperature perturbation as a function of depth 
and did not treat perturbations in the substrate below the layer. In contrast, our approach 
involves a computation of the entire eigenfunction from first principles, without any pre- 
conceptions as to the shape of the mode, and allows for perturbations in the substrate. It 
should be emphasized that we consider perturbations in all variables, including the compo- 
sition parameters X, Y and Z (—1 — X — Y). This enables us to find both thermal modes 
and nuclear modes (§6.2). In contrast, all previous studies restricted themselves to thermal 
perturbations, that too approximately, and thus were sensitive only to thermal modes. 

To get some sense of how good the approximate criteria of the past are, we present 
here some quantitative comparisons with two methods described in the literature, those of 
Fushiki & Lamb (1987a) and Cumming & Bildsten (2000). These papers use two different 
criteria for the onset of instability. Both express the criterion as 

~dl r> ^F' (38) 

but differ in their definitions of the quantities. Fushiki & Lamb (1987a) define eh ea t to be 
the total nuclear energy generation rate while Cumming & Bildsten (2000) define it to be 
1.9 times the energy generation rate from the triple-alpha reacton. In order to include the 
hydrogen bursts in the latter calculation we generalized it by adding also the hydrogen energy 
generation rate. The differences in the cooling rates in the two approaches are more striking. 
Fushiki & Lamb (1987a) asume that dTi/dH = and use a formula that is equivalent to 
(this is our best guess since their paper does not explain in sufficient detail) 

de coo i d ( dF\ d\wK ^ d 2 \nK ^ 



dT dT\dEJ dT dTdZ 

where e is the total energy generation rate and K = 16<rT 3 /3/t is the thermal conductivity 
(see eq. 4). Cumming & Bildsten (2000), on the other hand, start from the work of Fujimoto 
et al. (1981) and write approximately 



(acT A \ 
dT ~ dT \3kE 2 J ' 



d^cooi d 



(40) 
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where the quantity on the right is roughly the thermal energy density divided by the diffu- 
sion time to the surface. Clearly both criteria are approximate. Both are also local since 
they focus on values at a particular depth. The Fushiki & Lamb (1987a) method could be 
generalized to an approximate global criterion (it is possible that the authors did make such 
a generalization, but we have not been able to understand exactly how they might have done 
it). We limit ourselves to the local version given above. The Cumming & Bildsten (2000) 
criterion is of a mixed kind, since it compares a local quantity for the heating rate to a 
pseudo-global quantity for the cooling rate. 

We use the equilibria as computed with our code using our particular prescriptions for 
opacity, equation of state, etc. (§2.3) for all the calculations. We thereby ensure a fair 
comparison of the different methods. We consider the calculations done with our full linear 
perturbation analysis as the "correct" answer. Indeed, since the two approximate criteria 
described above are both restricted to thermal perturbations, we feel that the "correct" 
answer is the set of thermal modes calculated with our perturbation analysis (filled circles 
in Figs. 16 and 17). 

Figure 16(a) shows results corresponding to a temperature boundary condition, with the 
crosses depicting the results obtained with the Fushiki & Lamb (1987a) criterion. For the 
latter, we applied the local criterion (38), (39) at the temperature maximum in the accreted 
layer. (By choosing this point rather than the bottom of the layer, we feel that the method 
has a better chance of approximating the true global problem.) The agreement between the 
approximate criterion and our full calculations (thermal plus nuclear modes) is quite good. 
The Fushiki & Lamb (1987a) criterion finds a smaller range of / acc for the helium bursts and 
overestimates the range of the hydrogen bursts. Also, not surprisingly, it does not find the 
delayed mixed burst regime. The most surprising result is that the method, which is overtly 
limited to thermal perturbations, finds modes well outside the range of thermal modes. This 
may be viewed either as a serious error in the method or as an unexpected strength! 

Figure 17 shows the corresponding comparison for the Cumming & Bildsten (2000) cri- 
terion. Since these authors use a flux inner boundary condition, we calculated the equilibria 
as well as our stability results for this boundary condition, assuming an outward flux of 150 
keV per accreted baryon at the bottom of the accretion layer. In the stability analysis we 
assumed that the perturbed flux is zero at the bottom of the layer (which corresponds to a 
non-zero temperature perturbation). Once again, we should view the thermal modes (filled 
circles) as the "correct" answers. The crosses show the range of instability according to 
the criterion (38), (40) applied at the bottom of the accreted layer (where the temperature 
reaches a maximum). Wherever the Cumming & Bildsten (2000) criterion indicates that 
there is an instability, the agreement in the value of £i ay er,crit is fairly good. However, the 
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Fig. 17. — Models for log R/R s = 0.4 and M = 1.4M with a flux inner boundary condition. 
It is assumed that a constant flux of 150 keV per accreted baryon flows out of the star into 
the accretion layer. The open circles correspond to the full calculation where all variables are 
perturbed, and the filled circles to a calculation in which the composition is not allowed to 
vary (thermal modes). The crosses trace the results corresponding to the instability criterion 
of Cumming & Bildsten (2000). 
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method appears to miss finding an instability over a fairly wide range of / acc . Also, naturally 
the method does not reproduce the delayed mixed bursts (which we have identified as an 
effect associated with nuclear modes). 

One comment is in order. In calculating the above results for the Fushiki & Lamb 
(1987a) and Cumming & Bildsten (2000) criteria, we followed the philosophy described in 
§3, namely to label a system unstable only if it has unstable configurations on the "wall" of 
the S-curve for values of Ei ayer above the hydrogen and helium peaks. Thus, a situation like 
Fig. 2(a) in which the wall solutions are all stable is considered stable, even though there 
are unstable solutions on the trailing slope of the helium peak. We do not know if the other 
authors used this criterion or not, but we believe it is the correct approach to the problem 
(for reasons discussed in §3). 

In summary, the two approximate criteria discussed above are useful for quick estimates, 
but are not appropriate for detailed quantitative results. For the latter, one needs to do a full 
linear stability analysis, as presented in this paper, or carry out time-dependent simulations 
of the kind reviewed in §1. 

The good news is that, in the regions where both our and previous work have found 
instability, the results agree quantitatively quite well. Specifically, the results presented by 
Fujimoto et al. (1987b) agree in detail with the equilibrium curve presented in Fig. 4(a). For 
low accretion rates and cool cores, Hanawa & Fujimoto (1986) and Fujimoto et al. (1987a) 
found very long X-ray bursts with durations similar to those shown in Fig. 13. We also find 
good qualitative agreement with the more comprehensive results of Fushiki & Lamb (1987a) 
and the different burning regimes that they have identified, and quantitative agreement 
with the tabulated column densities and temperatures in Table 2 of Cumming & Bildsten 
(2000). However, the overarching nature of the calculations presented here make it difficult 
to compare very precisely with earlier work which were generally focused on understanding 
specific ignition criteria, rather than presenting predictions of the gross properties of Type I 
bursts over a wide range of accretion rates, temperatures and stellar radii. 



6.4. Comparison with Observations 

6.4-1 ■ Stable Burning Limit 

Van Paradijs et al. (1988) have collected together a body of burst data obtained with 
the EXOSAT satellite. They define a quantity 7 to be the ratio of the observed persistent 
flux F p to the net peak flux F re , not including the persistent emission, of bursts with radius 
expansion. The quantity F p is proportional to L acc , while F re is proportional to L' Edd , where 
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L' Edd is the Eddington luminosity in the frame of the neutron star surface and is smaller 
than our L Edd (defined at infinity) by a redshift factor. Thus, 7 is larger than L acc /L Edd by 
a variable factor that may be ~ 1.5. Van Paradijs et al. (1988) find that bursts occur for 
l°g 7 ^ —0.5; sources with 7 above this limit apparently burn fuel stably without bursting. 
The observational limit they have obtained on 7 should be compared with our prediction 
that bursts occur only below log (L acc /LEdd) ~ —0.7 to —0.6 (§4.3). The agreement is fairly 
good. In a previous paper (Narayan & Heyl 2002), we found that bursts should occur up 
to log(Z acc ) ~ —0.5. However, those calculations were done with a less sophisticated inner 
boundary condition (see §2.2). 

The brightest low-mass X-ray binaries generally do not show burst activity. Matsuba et 
al. (1995) saw Type I bursts from the bright source GX 13+1 when 7 ~ 0.3, Kuulkers & van 
der Klis (2000) saw a radius-expansion burst from the source GX 3+1 when the persistent 
luminosity was 0.17L E dd, and Kuulkers, van der Klis & van Paradijs (1995) and Smale (1998) 
observed bursts from Cyg X-2 with 7 as large as 0.74. Tournear et al. (2003) have recently 
used the USA satellite to carry out long-duration observations of a number of neutron star 
systems. They observed bursts from sources with a range of / acc extending up to ~ 0.3. 
Finally, Cornelisse et al. (2003) find that bursts cease at a luminosity of 5.5 x 10 37 ergs -1 , 
which corresponds again to 0.3L Edd (for M = 1.4M Q and our definition of the Eddington 
luminosity). These observations are all generally consistent with the van Paradijs et al. 
(1988) results. The stable limit according to our model is encouragingly close to the data, 
which is noteworthy since our model is a first-principles calculation with no free parameters. 

6.4.2. Burst a 

As mentioned in §5.4, both the predicted recurrence time and the fluence of bursts 
are subject to an uncertainty over whether the entire reservoir of fuel is consumed during 
the burst. Only fully time- dependent calculations can predict what fraction of the fuel is 
consumed in the burst, and even then probably only with multi-dimensional simulations. 
This is well beyond the scope of our calculation. However, partial burning affects the recur- 
rence time and the burst fluence in an identical fashion, and therefore their ratio should be 
insensitive to this uncertainty. The parameter a defined in §5.4 is such a quantity. 

Van Paradijs et al. (1988) present results for a as a function of 7. Fig. 18 compares 
their data with our predictions. In plotting the data, we have assumed that the observed 7 is 
equal to Z acc . As discussed in §6.4.1, the two may differ by a few tens of percent (because of 
gravitational redshift), so the data may need to be shifted to the left by this amount. Allowing 
for this uncertainty, we find that the data agree with the model predictions fairly well for 
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Fig. 18. — A comparison of the model predictions for a with the observational data presented 
in van Paradijs et al. (1988). The circles show the data. The lines in each panel correspond 
to the model results shown in Fig. 15, with the symbols replaced by connected lines. In each 
panel, the uppermost line corresponds to a neutron star radius of 6.5 km, the middle line to 
10.4 km, and the lowest line to 16.4 km. 
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luminosities above log(/ acc ) ~ —1.5. For accretion luminosities in the range log(/ acc ) ~ —1.5 
to —1, the data give a ~ 10 2 . This value is roughly consistent with our model predictions 
for smaller neutron star radii (6.5 km and 10.4 km); the largest radius we have tried (16.4 
km) does not agree very well with the data. 

The most interesting feature of the data is that, for higher accretion luminosities, a 
shoots up to values ~ 10 3 . Van Paradijs et al. (1988) suggested that the increase is because 
much of the nuclear energy in the accumulating fuel is released before the configuration 
becomes unstable, i.e. between bursts during accretion. The regime of increased a in the 
data very nicely coincides with the regime of delayed mixed bursts in our models. Indeed, 
in this regime, much of the nuclear fuel does burn steadily, as suggested by van Paradijs 
et al. (1988), and it is only after a considerable surface density has accreted that the burst 
is initiated. Correspondingly, a goes up by a large factor. According to our models, in this 
regime, the critical £i ayeriCri t to initiate a burst increases (Fig. 12), as does the recurrence time 
between bursts (Fig. 11), whereas the burst fluence as measured by tH+He is either unchanged 
or goes down (Fig. 13). All of these features are consistent with the observations, and this 
is one of the successes of the present work. 

The sudden "break" from values of a ~ 100 to larger values occurs at log(Z acc ) ~ — 1 
in the van Paradijs et al. (1988) data, though the precise location is hard to determine. 
A similar break has been seen in recent data presented by Cornelisse et al. (2003) who 
quote a break luminosity of 2.1 x 10 37 ergs -1 (they actually present data for burst rates 
and burst durations, but it is straightforward to translate these to a). In our units, their 
break occurs at log(/ acc ) = —0.9, which is close to where we find a break in our models. This 
quantitative agreement is very encouraging. Cornelisse et al. (2003) interpret the break in 
terms of a transition from unstable to stable hydrogen burning, even though such a transition 
is expected to occur at very low luminosities, not at log(/ acc ) ~ —1, and it does not predict 
the particular behavior seen in the data. In our model, the break signals the switch from 
the regime of prompt mixed bursts to that of delayed mixed bursts. Both the position of the 
break and the signatures we predict are in encouraging agreement with the observations. 

At low accretion luminosities ~ 0.01L Edd , van Paradijs et al. (1988) find very low values 
of a ~ 10. We are unable to reproduce this trend and are not aware of any other studies that 
succeed. Our models predict an increase in a as we move into the regime of helium bursts 
(where there is little or no hydrogen to burn), which is just the opposite of what is seen in 
the data. It would be useful to probe this regime in more detail with modern observations. 
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Fig. 19. — A comparison of the model burst durations tH+He with the corresponding data 
presented in van Paradijs et al. (1988). The model values are from Fig. 13, with the symbols 
replaced by connected lines. 
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6.4-3. Recurrence Time 

The results on t TCC in Fig. 11 show a variety of patterns that could be tested against 
observations, but we are not aware of appropriate data in the literature. Continuous moni- 
toring of sources, e.g., the recent work of Tournear et al. (2003) and Cornelisse et al. (2003), 
would be very useful in this regard. The latter paper finds that the peak burst rate in several 
sources is about 10 bursts per day, corresponding to a recurrence time of about 2.5 hours. 
(They also see sources with a peak burst rate of 2.5 per day, but these seem to belong to 
a separate class.) In our model, the shortest recurrence time is about 5 hours (Fig. 11), 
corresponding to a burst rate of about 5 bursts per day; this particular rate is found for 
smaller neutron star radii (6.5 km and 10.4 km) at the break point between the prompt 
mixed bursts and delayed mixed bursts at log(/ acc ) ~ —1. The location of the peak is in 
good agreement with the observations, but there is a factor of 2 discrepancy in the maximum 
burst rate. The latter may indicate that some of the approximations we have made (e.g., 
the opacity or the inner boundary condition) are still inadequate. It might also indicate that 
only a fraction of the fuel is burned during the burst. 

The observations of Cornelisse et al. (2003) indicate in a few sources that the burst rate 
decreases proportional to the luminosity below the break point (e.g., GX 354-0 in their Fig. 
2); equivalently, their recurrence times vary inversely as t rcc oc l~^ c . Such a trend corresponds 
to a constant value of S layer crit , and arises naturally in our model (and in other models, e.g., 
see Figs. 16, 17). Among the models we have calculated, it appears that smaller stellar radii 
(6.5 km, 10.4 km) are more consistent with the observed trend than larger radii (16.4 km). 
Also, of the four core temperatures we have tried, we find best agreement with the data for 
the hotter temperatures 10 7 ' 5 , 10 s , 10 8 5 K. A core temperature of 10 7 K appears unlikely. 

Above the break, Cornelisse et al. (2003) find a decrease in the burst rate by a factor 4. 
This is consistent with the increase in t rec seen in the models. Another interesting feature 
to check would be the sudden increase of t rcc with decreasing Z acc in the helium burst regime 
(^acc < —1-5). This regime has regions of very long recurrence times, which for all practical 
purposes may be considered burst-stable since observations are unlikely to detect these very 
rare bursts. We thus predict a "gap" in the occurrence of bursts for accretion luminosities 
in the vicinity of ~ 10~ 2 L E dd (Narayan & Heyl 2002). It is not that there are no bursts 
here, but that bursts are very rare and correspondingly very luminous (Fig. 13). It is not 
clear that the helium burst regime has been seen in the observations, or that there is a gap. 
Indeed, the van Paradijs et al. (1988) give anomalously low values of a here, whereas helium 
bursts should have large values of a. 
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6.4-4- Burst Fluence 

Figure 19 compares observed burst durations r from the van Paradijs et al. (1988) 
compilation with our model predictions for the Eddington-scaled duration tH+He- We see 
that the predicted durations are generally larger than the observed ones, though the models 
with a stellar radius of 6.5 km (the lowest lines in the panels) come pretty close. As in 
the case of the recurrence time, this comparison again is compromised by the possibility 
of partial burning of the fuel. Modulo this caveat, we conclude once again that somewhat 
hotter cores T core > 10 7 5 K and smaller neutron star radii R < 10 km are preferred. 

One interesting feature in the models is that the burst duration drops by about a factor 
of 1.5-2 above the break point at log(/ acc ) ~ —1. This is precisely the same point at which 
the burst recurrence time increases by a large factor. In the data presented by Cornelisse 
et al. (2003), such a decrease in the burst durations is seen (it is less obvious in the van 
Paradijs et al. data), but the amount of the decrease appears to be larger than we predict. 

6-4-5. Overstable Modes 

Revnivtsev et al. (2001) have reported a new class of low-frequency oscillations in burst 
sources that might be related to oscillatory behavior of the accreted layer prior to a burst. 
The oscillations seem to be present for accretion luminosities of about a tenth Eddington, 
and the frequency is about 0.01 Hz. 

Because our analysis computes the complex frequencies of unstable modes, it is natural 
to ask if we can explain the observed oscillations. The oscillations seem to be observed mostly 
in the regime that we have identified as delayed mixed bursts, above the break discussed in 
the previous subsections. In this regime, we do find complex values of 7 with the imaginary 
parts often larger than the real parts. However, our mode frequencies are generally lower 
than the observed frequencies by a factor of 10 or more. We do not know if this is an 
indication that some of our input physics is inadequate (§6.1) or perhaps that our numerical 
technique for finding eigenmodes is unable to converge on modes with very high frequencies 
(high compared to the accretion frequency 7 acc )- More work is needed. Alternatively, the 
oscillations may be related to slow propagation of the burning front on the stellar surface 
(see Bildsten 1995). 
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6.5. Future Prospects 

The richness of the properties of nuclear burning on neutron stars as revealed by the 
ab initio linear stability analysis presented here is remarkable. Considering the encouraging 
agreement that we find with observational data, the model deserves further exploration. 
Since we have used a variety of prescriptions for calculating the opacity, the equation of 
state, the nuclear burning, etc., it might be useful to improve these aspects of the model. 
The radiative free- free opacity in the outer layers seems to have a large effect (§6.1), so this 
is one area that needs more work. The treatment of the inner boundary condition is also a 
serious issue that needs to be considered in greater detail, especially for cool neutron star 
cores (see §2.2). 

One simplification in the present calculations is that we consider only a limited number 
of reactions for H-burning and He-burning (§2.3.3). At a sufficiently high temperature, 
breakout from the hot CNO cycle is expected to occur and hydrogen will be burned via the 
rapid proton capture (rp) process (Wallace & Woosley 1982; Schatz et al. 1999). Figure 
5 of Schatz et al. (1999) shows the temperature limits above which the reaction rates 
for processes such as 15 0(a, 7) 19 Ne and 14 0(a,p) 17 F exceed the usual /3-decay and proton 
capture rates for these nuclei. We have confirmed that all the burst-unstable configurations 
we have calculated in the present paper lie to the left of the corresponding lines in the Schatz 
et al. plot. Hence, none of our burst results are affected by the neglect of these additional 
reactions. However, close to the Eddington accretion rate, where our model predicts that the 
systems should be burst-stable, some of our equilibria cross into the zone where the breakout 
reactions, especially 15 0(o;, 7) 19 Ne, become important. It is thus conceivable that some of 
the very high M configurations that we have identified as stable may turn out to be unstable 
when breakout reactions are included. This topic needs further investigation. 

In terms of further calculations, the most obvious next step is to explore the effect 
of varying the composition of the accreting gas. Variations in metallicity between burst 
sources in globular clusters and those in the Galactic disk may possibly have observational 
consequences. As a more extreme example, ultracompact binaries with degenerate helium 
secondaries should have very different properties from the hydrogen- dominant systems we 
have considered here. Obviously, without hydrogen, there will be no unstable hydrogen 
burning, so the hydrogen burst region will be missing from the results. Instead, we imagine 
the helium burst regime will extend both to lower and higher accretion rates than in the 
results presented in this paper. Because helium requires a higher temperature/density to 
ignite, the bursts should occur more rarely than for hydrogen-accreting neutron stars. 

In Narayan & Heyl (2002), we examined the burst properties of more massive primaries 
(black hole candidates) and found that the threshold for stable nuclear burning moves to a 
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higher accretion rate. Hanawa & Fujimoto (1982) argue that the properties of the nuclear 
burning to lowest order are functions of the surface gravity and the mass accretion rate; 
therefore, one can understand the mass dependence of burning regimes by keeping the mass 
fixed and varying the radius and accretion rate; specifically, increasing the mass of the pri- 
mary shifts a diagram like Fig. 10 down and to the right which helps to account for the results 
of Narayan & Heyl (2002). Even though this argument gives a qualitative understanding of 
the mass dependence of the burning regimes, serious quantitative applications requires the 
kind of detailed calculations we have presented here. 

Recently, a class of very long X-ray bursts has been discovered (Cornelisse et al. 2000; 
Wijnands 2001) which have been interpreted as carbon-burning bursts (Taam & Picklum 
1978; Strohmayer 2000; Cumming & Bildsten 2001; Strohmayer & Brown 2002). In most 
cases, the bursts cannot be due to helium burning of a very thick layer since normal Type 
I bursts are seen at the same accretion rate, sometimes just before the superburst occurs. 
It would be interesting to include carbon-burning physics in our model to see what kinds of 
bursts are predicted. It is worth noting, however, that not all long bursts necessarily arise 
from carbon-burning. As we see in Fig. 13, we predict very long helium and even hydrogen 
bursts at low accretion luminosities under appropriate conditions. Some of the observed long 
bursts (e.g., Gotthelf & Kulkarni 1997) may correspond to these. 

On the observational front, there are numerous tests that one could envisage based on 
the results shown in Figs. 10-15. The different regimes of bursting — delayed mixed bursts, 
prompt mixed bursts, helium bursts, hydrogen bursts — reveal very distinct patterns in 
various observables. If these patterns are seen in the data, then one might be able to 
constrain the neutron star radius and/or the core temperature fairly well. This would have 
important implications for the neutron star equation of state and the nature of neutrino 
cooling in the core. We have made a beginning along these lines in this paper and have 
argued that perhaps T corc is > 10 7 5 K and neutron star radii are < 10 km. More work along 
these lines is worthwhile. 

If we are to extend this model, which has been developed for neutron stars, to ther- 
monuclear instabilities in accreting white dwarfs (classical novae), it would be necessary 
to include Coulomb corrections in the equation of state, a more accurate treatment of the 
semi-degenerate regime, and more detailed radiative and conductive opacities. Understand- 
ing whether there is a regime of stable nuclear burning at high accretion rates onto a white 
dwarf is a key ingredient of any scenario in which Type la supernovae result from stable ac- 
cretion of mass onto a white dwarf (Hillebrandt & Niemeyer 2000). Because classical novae 
typically eject the material accreted along with some of the substrate, a regime of stable 
nuclear burning is required for accretion to cause the white dwarf to grow in mass and to 
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end up in a supernova. Unfortunately, our current prescriptions are too crude to treat this 
important burning regime accurately, but the potential of the technique for white dwarfs is 
tantalizing. 

Finally, the methods that we have described here have deliberately avoided examination 
of the physics of the burst itself, during which many complications arise, including other 
nuclear reaction channels, convection and hydrodynamics. An important question to examine 
is whether our techniques could be extended in any simple way to study the properties of 
the bursts. 



7. Summary 

We have presented a comprehensive treatment of the stability of nuclear burning on 
the surface of neutron stars. For the first time, we have calculated the linearly unstable 
eigenmodes of an accretion layer in quasi-steady state, making no ad hoc assumptions re- 
garding the criterion for instability. The model reproduces the various previously known 
regimes of nuclear burning on neutron stars, and agrees with earlier results where there is 
overlap. Additionally, we have been able to probe in detail the behavior of accreting neutron 
stars at high mass accretion rates, near the threshold of stable nuclear burning. Here, we 
find a hitherto unrecognized regime of delayed mixed bursts, with very distinct properties 
compared to the more standard prompt mixed bursts. 

For accretion rates greater than one percent of the Eddington rate, we find encouraging 
agreement between the model predictions and observations of bursts. The existence of the 
regime of delayed mixed bursts provides a natural explanation for the observed dramatic 
increase of the burst parameter a at high accretion luminosities (Fig. 18 and Cornelisse et 
al. 2003). In addition, there is some indication from the preliminary comparisons presented 
here that burst systems have hot cores > 10 7 ' 5 K, consistent with cooling in the neutron star 
interior being dominated by the modified URCA process or a similar low-efficiency cooling 
mechanism. Cool cores with T core ~ 10 7 K, as might be present if direct URCA cooling were 
to operate, are less likely. We also find a number of indications for small neutron star radii 
< 10 km. These results could be tightened with more careful modeling, e.g., by improving 
some of the prescriptions we use for the input physics (§6.5), and with more extensive and 
better quality data. 

We thank Deepto Chakrabarty, Duncan Galloway, Erik Kuulkers, Feryal Ozel and Dim- 
itrios Psaltis for useful discussions and comments. We are also grateful to the referee for a 
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Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for 
and on behalf of NASA under contract NAS8-39073. 



A. Equilibrium and Perturbations 

Here we discuss in more detail the basic equations describing the equilibrium and per- 
turbations of the accreted layer. The governing equations are given in equations (3)-(7). 
We use the notation defined in Table 1 and in equations (14), (15), where quantities like 
Po, Tq refer to the equilibrium, p, T refer to the corresponding quantities in the perturbed 
state, and pi, T\ refer to the spatial component of the linear perturbations. The equilibrium 
quantities have no time dependence, while the perturbations have a time dependence of the 
form exp(7t). Since 7 is in general complex, the quantities p, T, p 1 , 7\, etc. are complex, 
whereas p , T , etc. are real. Also, since X + Y + Z = 1 and X + Y + Z = 1, we have 
X 1 + Y 1 + Z 1 = 0. 

Let us begin with the radiative transfer equation (4). For the equilibrium, this equation 
gives 

dT _ 3k F 

<9S 16aT 3 ' 1 ' 

where all quantities are real. For the perturbations, we linearize equation (4) and consider 
first-order deviations. This gives 



<9S 16<7T 3 1 16a 



dp 
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/ 
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where now the various quantities are in general complex. For instance, d(n/T 3 )/dp refers 
to the derivative of the complex quantity k/T 3 with respect to complex p (the derivative is 
well-defined since all the quantities are analytic). In practice, we calculate such derivatives 
numerically. Since k in general depends on all three quantities X, Y, Z, it is necessary to 
replace Z by 1 — X — Y or Z\ by —X 1 — Y\ before computing the partial derivatives d/dX, 
d/dY. 

Instead of writing an equation for the linear perturbation 7\, we could equally well 
consider the equation for the total (complex) perturbed temperature T = T + T 1 exp(^t)Ti. 
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This is nothing but the original equation 

dT 3kF 



as 16crT 3 ' 



(A3) 



where F is now the total complex perturbed flux, and k is the opacity corresponding to 
the perturbed p, T, X, Y. Within the linear approximation, this equation, coupled with 
equation (Al), has the same content as equation (A2) for the linear perturbations (equation 
A2 multiplied by exp(7t) is just the difference of equations A3 and Al). Equation (A3) 
has the advantage of being more compact than equation (A2). It is also numerically more 
convenient, since the compactness of the equation translates to relative simplicity of the 
corresponding computer code. 

Consider next the H-evolution equation (6), which involves a time derivative. The 
equilibrium is described by 

• 8X en,o 



<9E 



where all quantities are real. The linear first order perturbation X\ satisfies 



,x 1 + t 9x ^ 
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where the term 7X1 on the left comes from the time derivative d/dt operating on exp(7t)Xi. 
Once again, instead of considering the equation for the perturbation Xi, we may write down 
the equation for the total perturbed quantity X = X + exp(-yt)Xi. 

7( X_X 0) + £_ = __|, (A6) 

which is nearly identical in form to equation (A4) for the equilibrium, except that (i) it has 
an extra term proportional to 7 because of the time dependence of the perturbations, (i) en- 
is evaluated at the perturbed p, T, X, Y, and (iii) all quantities are complex. 

Following the above examples, the other three equations are straightforward. The hy- 
drostatic equilibrium equation (3) gives for the equilibrium 

dPpdpo dPpdTo dP dX dP dY = 

d Po as dT as ax as ar as 9: [ ' 

and for the perturbations 

dP dp dPdT dPdX dP dY _ 
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The energy conservation equation gives 



dF 
<9£ 



= -ST n 



— — = -7T 
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Finally, the He-evolution equation gives 

^ dY = eufi 
<9£ ~ £* 

• dY 

Equations (Al), (A4), (A7), (A9) and (All) are five differential equations for the equi- 
librium quantities F , p , T , X , Y . We assume a value for the outgoing flux F out at the 
surface of the layer, solve for p ou t,o, T ontt0 , X outi o, ^ O ut,o from the outer boundary conditions, 
and then integrate the 5 differential equations down to the bottom of the layer and then into 
the substrate to a depth equal to S diff defined in equation (18). At this depth, we require 
the temperature T (S diff ) to be equal to the required core temperature T core . If it is not, 
we change the value of -F O ut,0 and repeat until the inner boundary condition is satisfied. We 
then have a valid equilibrium solution. 

For the perturbations, we work with equations (A3), (A6), (A8, (A10), (A12), which are 
five differential equations for the total perturbed quantities F, p, T, X, Y. At the surface, 
we set the flux equal to F out = F outi o + i^out.i, where F outj i <C F outt0 in order to satisfy the 
assumption of a linear perturbation of the equilibrium. We solve for the other variables at 
the surface, assume a value for the eigenvalue 7, and integrate the equations down to a depth 
S max in the substrate (see eq. 20). At this depth we require T(S max ) = T (E max ). We vary 
7 until this inner boundary condition is satisfied, at which point we have a solution for the 
complex eigenvalue 7, and also the shape of the eigenfunction (iq, pi, etc.). The search 
in 7 space is tailored to find eigenvalues with positive real parts since only such modes are 
unstable. 

Let us define a turning point in the sequence of equilibria (S-curve) to be a point at 
which the locus of equilibrium solutions satisfies the condition (Ei ayer /dF out) o = 0. Since the 
derivative is zero, it means that two equilibria with escaping surface fluxes equal to F outj0 and 
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F ou t = F outj0 + F outj i, where F outj i <C F outt0 , both satisfy the differential equations as well as 
the boundary conditions for the same value of Si ayer . Now, the solution for F out satisfies the 
equilibrium equations (Al), (A4), (A7), (A9), (All), while the solution for F out satisfies the 
perturbed equations (A3), (A6), (A8), (A10), (A12). The only difference between the two 
sets of equations is the presence of various terms involving 7. If the second solution is also 
an equilibrium solution, then it implies that it must satisfy the perturbation equations with 
7 precisely equal to 0, i.e., one of the modes of the system corresponds to frequency. This 
result is not surprising. Since the initial model is at a turning point, there are neighboring 
equilibria with the same value of £i aycr but different values of F out; o, i.e., the system has a 
linear zero-frequency mode. 

The above theorem, that a mode with 7 = exists at turning points of the S-curve, 
is well-known for one-zone models, e.g., see the discussion of Paczynski (1983a) for an ap- 
plication to bursts. Our discussion shows that the same result is valid even when one is 
considering the more complex model described here. The reason it works is that our equi- 
libria are ultimately labeled by a single parameter, the value of -F ut,o at the surface. This 
is all that matters, and the fact that our solutions involve many variables and are described 
by continuous functions is not relevant. 
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